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Abstract 


Objectives  Soot  formation  in  aircraft  engines  is  of  great  concern  due  to  the  health  and 
environmental  impacts  of  soot.  The  design  of  aircraft  propulsion  systems  would  greatly 
benefit  from  a  computational  methodology  that  can  accurately  predict  soot  emissions.  This 
project  has  sought  to  develop  these  simulation  capabilities  in  the  framework  of  large  eddy 
simulation  (LES).  The  complexity  of  the  soot  formation  process  and  its  strong  interactions 
with  gas-phase  combustion  chemistry  and  turbulence  necessitated  a  comprehensive  approach 
involving  three  key  research  areas:  chemical  modeling  of  soot  surface  reactions,  especially 
oxidation,  statistical  modeling  of  soot  particle  distributions,  and  soot  modeling  in  LES  of 
turbulent  combustion. 

Technical  Approach  Soot  growth  and  oxidation  processes  were  investigated  using  theo¬ 
retical  approaches  applicable  to  large  reaction  systems.  The  work  focused  on  the  thermody¬ 
namics  and  kinetics  of  aromatic  oxyradicals,  the  key  intermediates  in  oxidation  of  soot. 

A  novel  moment  method  called  the  hybrid  method  of  moments  (HMOM)  was  developed 
for  predicting  moments  of  the  soot  particle  number  density  function.  This  formulation 
described  soot  particle  size  using  both  volume  and  surface  area  and  is  able  to  account  for 
bimodal  distributions  of  soot  particles  sizes. 

An  integrated  LES  model  for  soot  evolution  in  turbulent  nonpremixed  flames  was  de¬ 
veloped.  The  model  components  include  a  soot  model  based  on  HMOM,  an  extended 
flamelet/progress  variable  gas-phase  combustion  model,  and  a  model  to  account  for  the  slow 
chemistry  of  soot  precursors.  The  development  of  these  modeling  approaches  was  aided  by 
the  results  of  a  first-of-its-kind  direct  numerical  simulation  (DNS)  study  of  soot  nucleation 
and  growth  using  a  polycyclic  aromatic  hydrocarbon  (PAH)  based  model  for  soot  inception. 
Results  The  principal  conclusion  of  the  soot  oxidation  study  is  that  only  the  outer  rings  of 
aromatic  oxyradicals  are  able  to  undergo  oxidation,  whereas  the  inner  rings  resist  oxidation 
in  combustion  environments.  We  also  found  that  the  variability  of  both  thermodynamic 
stability  and  decomposition  rate  of  aromatic  oxyradicals  can  be  explained  and  correlated 
with  substrate  aromaticity.  These  findings  will  support  the  development  of  practical  rules 
for  predicting  oxidation  rate  coefficients  for  aromatics,  thereby  enhancing  models  of  soot 
oxidation. 

The  LES  soot  model  was  validated  in  a  laboratory  scale  natural  gas  jet  diffusion  flame. 
Compared  to  experimental  measurements,  the  LES  provided  a  reasonable  prediction  of  the 
maximum  soot  volume  fraction.  Factors  extrinsic  to  the  soot  model  formulation  were  found  to 
be  major  sources  of  error  and  uncertainty  in  soot  prediction.  These  factors  include  modeling 
of  small  scale  gas-phase  mixing  rates  and  specification  of  the  kinetic  mechanism  for  PAH 
formation  and  growth. 

Benefits  The  overall  outcome  of  this  project  is  a  significant  advancement  in  the  sophistica¬ 
tion  of  soot  modeling  approaches  for  LES.  Using  these  approaches,  the  superiority  of  LES 
to  Reynolds- averaged  Navier  Stokes  (RANS)  methods  for  soot  modeling  was  demonstrated. 
The  LES  model  permits  the  effects  of  turbulence  and  gas-phase  chemistry  on  soot  forma¬ 
tion  to  be  assessed  qualitatively.  Sources  of  error  in  soot  prediction  were  clearly  identified, 
allowing  future  modeling  studies  to  be  targeted  to  these  areas. 
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1  Objectives 

Soot  formation  in  aircraft  engines  is  a  major  concern  for  all  aircraft  engine  manufactur¬ 
ers.  While  current  emission  requirements  are  based  only  on  smoke  visibility,  quantitative 
requirements  on  soot  emissions  from  aircraft  engines  will  be  established  soon  due  to  the 
many  negative  impacts  of  soot  on  human  health  and  the  environment.  These  regulations 
will  not  only  stipulate  the  total  soot  emission  levels  but  other  soot  characteristics  as  well. 
For  instance,  it  has  been  realized  that  small  particles  are  much  more  carcinogenic  than  larger 
particles  indicating  that  the  size  distribution  is  a  critical  factor.  Soot  particles  also  act  as 
condensation  nuclei  for  ice  crystals  forming  cirrus  clouds,  which  influence  the  earth’s  climate 
through  radiative  forcing  effects.  The  cloud  nucleation  rate  is  linked  to  both  the  number  of 
soot  particles  and  their  structure.  From  these  examples,  it  is  evident  that  future  regulations 
will  specify  detailed  characteristics  of  permissible  soot  emissions  including  both  the  allowable 
size  distribution  and  particle  number  density. 

Novel  techniques,  such  as  lean  direct  injection,  have  been  proposed  to  reduce  soot  emis¬ 
sions  from  aircraft  engines.  Although  it  is  known  that  such  techniques  might  reduce  the  soot 
volume  fraction,  little  is  known  about  the  effect  on  particle  size  distribution  and  the  particle 
number  density. 

Accurate  modeling  of  soot  formation  is  very  difficult  due  to  the  complex  underlying 
chemical  and  physical  processes,  including  a  large  sequence  of  gas  phase  reactions  form¬ 
ing  polycyclic  aromatic  hydrocarbons  (PAH)  followed  by  particle  inception,  particle/particle 
interactions,  condensation  of  PAH  on  the  particle  surface,  and  soot  particle  growth  and 
oxidation  by  heterogeneous  reaction  with  different  gas  phase  species.  In  turbulent  flows, 
moreover,  we  have  shown  in  the  past  that  the  very  subtle  interaction  of  turbulent  transport 
and  molecular  differential  diffusion  effects  can  strongly  influence  predicted  soot  concentra¬ 
tions.  Although  the  effects  of  turbulence  are  very  important  for  the  evolution  of  particulates, 
limited  prior  work  has  been  done  in  this  area.  Several  empirical  or  semi-empirical  models  for 
soot  formation  have  been  proposed  in  the  past,  but  they  are  far  from  being  predictive.  The 
development  of  more  predictive  models  for  soot  emissions  and  soot  particle  size  distribution 
will  enable  the  computational  optimization  of  present  engine  designs  and  the  development 
of  new  low-emissions  combustion  concepts. 

The  major  objective  of  this  project  has  been  the  formulation  of  a  comprehensive  simu¬ 
lation  methodology  to  enable  accurate  and  detailed  prediction  of  soot  formation  in  aircraft 
engines,  especially  military-type  aircraft  engines.  These  simulation  methods  will  provide  a 
valuable  tool  for  engineering  design  and  the  enable  the  development  of  strategies  to  mitigate 
soot  formation. 

Large  eddy  simulation  (LES)  is  the  most  advanced  approach  available  for  simulating 
turbulent  reacting  fluid  flows  of  the  scale  and  complexity  of  the  flow  with  in  an  aircraft 
combustor.  Three  main  issues  were  identified  in  the  development  of  predictively  accurate 
LES  soot  modeling  capabilities. 

First,  the  soot  particle  population  evolves  under  the  effects  of  several  complex  physi¬ 
cal  and  chemical  processes,  including  particle  nucleation,  particle  growth  through  particle- 
particle  collisions  and  surface  reactions,  and  particle  erosion  through  oxidative  reactions. 
While  the  underlying  phenomena  of  the  nucleation  and  growth  processes  are  comparatively 
well  established,  soot  oxidation  process  are  less  completely  understood.  Fundamental  chem- 
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ical  modeling  studies  of  soot  oxidation  were  carried  out  under  the  project.  The  results  of 
these  studies  can  be  used  to  create  simplified  representations  of  the  oxidation  process  for 
inclusion  in  LES  computations. 

Second,  a  statistical  description  of  the  soot  population  must  be  developed.  Mathemat¬ 
ically,  this  description  is  available  in  the  form  of  the  soot  particle  number  density  function 
(NDF).  However,  the  NDF  must  be  approximated  for  compatibility  with  the  LES  numerical 
discretization.  Moment  methods  are  commonly  used  to  approximate  such  density  functions. 
Due  to  the  characteristically  bimodal  distribution  of  soot  particle  sizes,  special  care  is  re¬ 
quired  in  applying  moment  methods  to  the  soot  NDF  to  ensure  that  the  resulting  statistical 
description  is  numerically  robust  and  physically  accurate.  A  novel  moment  method  was 
developed  under  the  program  to  meet  these  requirements. 

Third,  the  small  scale  interactions  of  soot,  turbulence,  and  gas-phase  chemistry  are  un¬ 
resolved  in  LES  and  must  be  modeled.  It  must  be  recognized  that  soot  formation  inside 
turbulent  aircraft  combustors  is  a  highly  intermittent  and  unsteady  process.  In  practical 
combustor  designs  such  as  the  rich-quench-lean  (RQL)  operating  mode,  soot  is  generated  in 
fuel-rich  high  temperature  regions  and  is  subsequently  oxidized  after  dilution  with  air.  The 
amount  of  soot  actually  emitted  is  the  result  of  the  differences  between  these  generation  and 
oxidation  processes.  More  importantly,  the  generation  and  oxidation  processes  occur  with 
very  high  intermittency,  indicating  that  the  temporal  history  of  soot  precursors,  turbulence- 
related  strain,  and  combustion  processes  are  critical  in  determining  emission  levels.  All  of 
these  processes  are  highly  complex  and  interact  nonlinearly.  Therefore,  accurate  modeling 
of  unresolved  soot/flame/turbulence  interaction  is  essential  to  achieving  quantitative  pre¬ 
dictions  of  soot  formation.  Under  this  project,  new  modeling  methods  for  capturing  these 
interactions  were  developed  and  validated. 

Under  the  project,  significant  advances  were  made  in  each  of  the  three  areas  of  (i)  chemical 
modeling,  (ii)  statistical  modeling  ,  and  (iii)  simulation  and  modeling  of  soot /flame/ turbulence 
interaction.  The  remainder  of  this  Report  is  organized  as  follows:  Section  2  provides  the 
relevant  background  information  and  characterizes  the  outstanding  issues  for  accurate  simu¬ 
lation  of  soot  formation.  Section  3  describes  the  computational  methods  used  to  investigate 
these  issues  and  reports  on  the  new  models  developed.  Section  4  discusses  the  results  of 
the  simulations  performed  and  discusses  the  significance  of  their  findings.  Finally,  Section  5 
summarizes  the  achievements  of  the  research  and  indicates  areas  for  ongoing  work. 
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2  Background 

In  this  section,  we  will  first  provide  an  overview  of  the  present  understanding  of  soot  forma¬ 
tion  and  evolution  in  relation  to  our  own  work  in  this  field. 

2.1  Oxidation  of  Soot 

Chemical  evolution  of  graphene  edges  is  one  of  the  key  processes  in  the  formation,  growth, 
and  oxidation  of  soot  and  its  aromatic  precursors  [4-6].  Assuming  that  the  soot  particle 
surface  is  comprised  of  molecular  aromatic  sites  (i.e.,  graphene  edges)  and  invoking  chemi¬ 
cal  similarity,  it  was  postulated  that  soot  surface  reactions  can  be  thought  of  as  chemically 
analogous  reactions  of  polycyclic  aromatic  hydrocarbons  (PAHs)  [7,  8].  Growth  reactions, 
those  increasing  PAH  size  by  reactions  with  gaseous  species,  primarily  acetylene,  received 
immediate  and  continuing  attention  (see,  e.g.,  [4-6]  and  references  therein).  However,  it 
was  also  recognized  that  the  chemical  analogy  was  not  sufficient  to  describe  reactions  taking 
place  at  the  surface  and  steric  effects,  neighboring  sites,  and  substrate  size  must  also  be  con¬ 
sidered  [9].  The  latter  realization  increased  substantially  the  number  of  possible  elementary 
reactions;  for  instance,  the  latest  detailed  graphene-edge  growth  mechanism  is  composed  of 
42  elementary  reactions  [10]. 

The  oxidation  of  PAHs  and  soot  surfaces  has  received  lesser  attention.  The  initial  detailed 
models  of  soot  oxidation  invoked  two  principal  steps:  an  elementary  reaction  of  a  surface 
radical  with  O2  and  an  atomistically- unresolved  attack  of  OH  on  a  generic  surface  site  [8,  11]. 
Both  reactions  were  assumed  to  remove  one  C  atom  per  O  atom  of  the  gaseous  reactant  and 
both  assumed  to  form  a  “pristine”  surface  site  (Csurface-H  or  Csurface-)  as  the  substrate 
product  of  oxidation.  The  former  reaction  was  modeled  following  the  postulate  of  chemical 
similarity  by  phenyl  +  O2.  The  oxidation  by  OH  was  described  using  the  collision  efficiency 
determined  in  flame  studies  [12],  This  oxidation  model  has  been  widely  adopted  in  modeling 
studies,  recently  with  different  values  of  the  rate  coefficients  [13]. 

Considering  the  knowledge  gained  with  growth  reactions,  where  the  reaction  chemistry 
of  larger  aromatic  edges  turned  out  to  be  much  richer  than  that  of  the  analogous  small 
aromatic  species,  there  is  clearly  a  need  for  similar  exploration  of  possible  pathways  for 
oxidation.  This  was  set  to  be  the  objective  of  the  present  work.  We  began  this  process  with 
the  examination  of  the  thermodynamic  stability  of  the  key  intermediates,  graphene-edge 
oxyradicals,  located  at  both  zigzag  and  armchair  sites,  and  then  proceeded  with  the  kinetics 
of  graphene  oxyradical  decomposition. 

2.2  Statistical  Soot  Modeling 

Modeling  of  the  statistics  of  the  soot  particle  population  is  critical  for  the  accurate  descrip¬ 
tion  of  soot  evolution  and  prediction  of  the  particle  size  distribution  in  the  exhaust  of  an 
engine.  Two  prominent  characteristics  of  typical  soot  particle  populations  must  be  consid¬ 
ered  to  develop  an  appropriate  statistical  description  of  the  soot  number  density  function 
(NDF).  First,  the  size  distribution  of  soot  particles  is  frequently  observed  to  be  bimodal  in 
experimental  studies  [14,  15].  The  first  mode  consists  of  small  particles  formed  by  persistent 
nucleation  from  PAH,  and  the  second  mode  consists  of  larger  particles  which  have  grown  by 
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collision  and  surface  reactions.  Second,  while  particles  in  the  first  mode  are  overwhelmingly 
spherical  in  form,  the  shape  of  particles  in  the  second  mode  is  highly  variable,  including 
both  spherical  particles  [14,  15]  and  fractal  aggregates  [16],  which  are  particles  composed  of 
many  small  spheres  called  primary  particles.  These  complexities  in  shape  generally  require 
the  NDF  to  be  parametrized  by  multiple  quantities  such  as  volume  and  surface  area  [17-19] 
to  obtain  an  accurate  description  of  the  particle  population. 

The  evolution  of  this  complex  bimodal  NDF  is  governed  by  the  Population  Balance  Equa¬ 
tion  (PBE).  However,  the  high  dimensionality  of  the  PBE  (including  spatial,  temporal,  and 
two  or  more  particle  size  coordinates)  makes  its  direction  solution  intractable.  Therefore, 
alternative  methods  are  used  to  obtain  an  approximate  description  of  the  soot  particle  pop¬ 
ulation.  The  available  methods  can  be  divided  into  two  classes:  those  that  solve  for  an 
approximate  form  of  the  NDF  and  those  that  solve  for  a  set  of  mean  quantities. 

There  are  two  basic  methods  for  solving  for  an  approximate  form  of  the  NDF.  The  most 
accurate  statistical  modeling  method  is  Monte  Carlo  (MC)  simulation.  In  these  simulations, 
stochastic  processes  govern  the  evolution  of  a  representative  population.  The  benefit  of  MC 
simulations  is  a  very  accurate  prediction  of  the  complete  NDF.  However,  this  accuracy  comes 
at  a  cost,  and  the  application  of  MC  simulations  is  limited  to  simple  configurations  such  as 
homogeneous  reactors  and  one-dimensional  laminar  premixed  flames.  A  second  approach 
is  to  use  a  sectional  method  (SM).  For  these  methods,  the  NDF  is  discretized  in  bins.  An 
equation  for  the  number  of  particles  in  the  bin  is  solved  for  each  bin.  Like  MC  simulations, 
accurate  predictions  of  the  complete  (discretized)  NDF  can  be  obtained  but  at  significant 
computational  cost  since  many  bins  are  required  for  each  NDF  coordinate  to  obtain  good 
accuracy.  Unlike  MC  simulations,  sectional  methods  are  deterministic  and  can  be  explicitly 
coupled  to  gas-phase  chemistry. 

The  second  approach  to  the  statistical  modeling  of  soot  is  generally  described  as  one 
of  several  moment  methods.  In  these  methods,  equations  are  solved  that  yield  information 
about  mean  quantities  rather  than  the  NDF  itself.  Moment  methods  are  superior  in  terms  of 
computational  cost  compared  to  SM  and  MC  simulations,  but  complete  information  about 
the  population  is  lost.  An  additional  shortcoming  is  that  the  equations  for  the  moments 
are  unclosed.  Moments  that  are  not  directly  solved  for  are  required  to  compute  the  source 
terms  for  the  moments  that  are  solved  for,  thus  requiring  further  modeling.  At  the  outset 
of  the  project  two  options  in  moment  methods,  each  with  a  different  closure  model,  were 
available:  the  Method  of  Moments  with  Interpolative  Closure  (MOMIC)  [20,  21]  and  the 
Direct  Quadrature  Method  of  Moments  (DQMOM)  [17]. 

In  MOMIC  [20,  21],  the  equations  for  the  moments  are  closed  by  logarithmic  polynomial 
interpolation.  The  widely  used  simplified  version  of  MOMIC,  which  interpolates  based  on 
positive  moments  only,  was  found  to  be  unable  to  represent  the  bimodal  nature  of  the  soot 
NDF  under  some  conditions [18].  However,  more  comprehensive  formulations  of  MOMIC  [20], 
based  on  positive  and  negative  moments,  should  offer  better  accuracy  in  these  circumstances 
[21]- 

In  the  original  Quadrature  Method  of  Moments  (QMOM)  [22],  the  NDF  is  approximated 
by  a  series  of  delta  functions,  and  the  moments  are  approximated  by  Gauss  quadrature. 
DQMOM  [17]  is  the  generalization  of  QMOM  to  multivariate  distributions.  In  DQMOM, 
transport  equations  are  solved  directly  for  the  weights  and  locations  of  the  delta  functions. 
The  inherent  multi-modality  of  DQMOM  is  well  suited  to  capturing  the  persistent  nucleation 
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mode  [18].  However,  in  order  to  obtain  the  source  terms  for  the  transport  equations  for  the 
weights  and  locations  of  the  delta  functions  from  the  source  terms  of  the  moment  transport 
equations,  a  linear  system  must  be  inverted.  Depending  on  the  shape  of  the  NDF  at  any 
given  point  and  the  choice  of  moments  solved  for,  this  inversion  may  be  ill-conditioned. 

In  light  of  these  findings,  a  new  moment  method  was  developed  in  this  project.  The 
Hybrid  Method  of  Moments  (HMOM)  [23]  combines  the  inherent  multi-modality  of  DQMOM 
with  the  numerical  ease  of  MOMIC.  HMOM  uses  the  approach  of  DQMOM  for  the  first  mode 
(small  spherical  particles)  and  the  approach  of  MOMIC  for  the  second  mode  (large  fractal 
aggregates)  in  order  to  capture  the  influence  of  the  persistent  nucleation  mode.  Details  of 
the  formulation  of  HMOM  are  provided  in  Sec.  3.2.1.  HMOM  has  been  successfully  used 
to  predict  soot  formation  in  a  variety  of  laminar  and  turbulent,  premixed  and  nonpremixed 
flames  in  LES  and  DNS  [23-25]. 

2.3  Simulation  and  Modeling  of  Soot /Flame/Turbulence  Interac¬ 
tion 

The  small-scale  interactions  between  turbulence,  gas-phase  chemistry,  and  soot  have  a  pro¬ 
found  effect  on  the  soot  formation,  growth,  and  destruction  processes  in  turbulent  reacting 
flows.  Turbulence  affects  the  soot  field  in  two  ways.  First,  due  to  a  very  large  Schmidt 
number,  soot  is  confined  to  very  thin  structures  which  are  subsequently  stretched  into  long 
filaments  by  the  turbulent  eddies.  Second,  soot  is  formed  from  Polycyclic  Aromatic  Hydro¬ 
carbons  (PAH),  which  are  formed  only  in  regions  of  low  scalar  dissipation  rate.  The  result 
of  these  properties  is  a  both  spatially  and  temporally  intermittent  soot  field. 

The  use  of  large  eddy  simulation  (LES)  is  emerging  as  a  standard  practice  for  turbulent 
combustion  in  practical  devices  [26].  Since  LES  resolves  only  large-scale  features  of  the 
flow,  combustion  and  soot  formation  which  occur  exclusively  at  the  small  scales  must  be 
modeled.  In  particular,  the  small-scale  correlations  between  gas-phase  species  and  soot 
precursor  evolution  need  to  be  described. 

The  first  three-dimensional  LES  of  soot  evolution  in  nonpremixed  turbulent  combustion 
was  performed  by  El-Asrag  and  Menon  [27].  Their  modeling  approach  was  based  on  the 
inclusion  of  a  moment  method  for  soot  within  the  Linear  Eddy  Model  (LEM)  for  the  subfilter 
combustion  processes.  The  model  was  shown  to  reasonably  predict  the  soot  volume  fraction 
in  an  ethylene  jet  flame.  However,  while  successful  in  the  case  considered,  the  modeling 
strategy  relied  on  a  simplistic  soot  model  and  a  semi-empirical  acetylene  based  soot  inception 
model. 

Recently,  the  small-scale  interactions  between  soot,  turbulence,  and  chemistry  have  been 
investigated  using  DNS.  The  works  of  Lignell  and  coworkers  [28,  29]  highlighted  the  role  of 
transport  and  found  that,  while  soot  had  roughly  the  same  probability  of  moving  towards  or 
away  from  the  flame,  the  soot  concentration  was  highest  when  soot  moved  toward  the  flame. 
However,  a  semi-empirical  acetylene  based  soot  inception  model  was  used  in  their  work. 
More  recently,  Bisetti  et  al.  [24]  performed  the  first  DNS  of  soot  evolution  with  a  PAH  based 
soot  inception  model  under  this  project.  Results  from  this  study  are  provided  in  Sec.  4.2, 
but  some  key  findings  are  summarized  here.  In  contrast  to  previous  DNS  studies,  Bisetti 
et  al.  found  that  soot  concentrations  were  highest  when  soot  moved  away  from  the  flame. 
The  data  showed  that  soot  moving  toward  the  flame  was  oxidized  quickly  before  appreciable 
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soot  growth  could  occur,  while  soot  moving  away  from  the  flame  had  long  residences  time 
during  which  significant  growth  by  PAH  condensation  occurred.  In  addition,  PAH  was  found 
not  only  to  be  very  sensitive  the  scalar  dissipation  rate,  which  is  partially  responsible  for 
soot  intermittency,  but  also  to  exhibit  significant  unsteady  kinetic  effects  due  to  relatively 
slow  chemistry.  This  work  highlighted  not  only  the  important  role  of  PAH  in  turbulent 
nonpremixed  combustion  but  also  its  difficulty  in  modeling. 

As  discussed  in  Sec.  2.2,  the  soot  population  is  described  using  a  set  of  moments  of 
the  underlying  number  density  function.  In  this  context,  small-scale  correlations  between 
soot,  turbulence,  and  gas-phase  chemistry  are  best  represented  by  a  one-point  one-time  joint 
probability  density  function  (PDF)  of  gas-phase  species  and  soot  moments.  Options  for 
modeling  this  joint  PDF  fall  into  two  basic  categories,  presumed  PDF  modeling  approaches 
and  transported  PDF  modeling  approaches.  Computational  details  of  these  approaches  are 
provided  in  Sec.  3.2.3. 

Due  to  the  inherent  challenges  in  describing  soot  precursor  chemistry,  turbulence,  and 
combustion  processes,  not  much  focus  has  been  given  to  the  modeling  of  this  joint  PDF.  In 
LES,  the  subfilter  correlations  are  often  neglected.  In  this  case,  the  PDF  could  be  written 
as  the  product  of  the  marginal  PDFs  of  gas  phase  scalars  and  soot  moments.  However, 
given  the  high  sensitivity  to  soot  location  with  respect  to  the  flame,  this  approximation  is 
bound  to  introduce  large  errors.  Recently,  Mueller  [30]  have  developed  a  presumed  PDF 
approach,  where  the  independence  of  the  variables  is  still  assumed  but  the  marginal  PDFs 
of  the  moments  are  modeled  using  an  intermittency  related  term.  In  the  Reynolds-averaged 
Navier  Stokes  (RANS)  formulation,  Lindstedt  et  al.  [31]  have  directly  solved  the  joint- 
PDF  using  a  high-dimensional  transport  equation.  It  was  found  that  accounting  for  the 
correlations  between  gas-phase  scalars  and  soot  moments  introduced  large  differences  in  the 
soot  profiles.  Chandy  et  al.  [32]  have  used  the  PDF  approach  in  the  context  of  LES,  but 
with  a  simplified  model  for  the  soot  chemistry. 
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3  Methods  and  Models 


This  section  presents  the  methods  used  to  study  soot  growth  and  oxidation  processes  (Sec.  3.1) 
and  the  models  developed  for  simulation  of  soot  evolution  in  LES  (Sec.  3.2).  The  LES  mod¬ 
eling  discussion  focuses  on  the  novel  closures  developed  in  this  project,  which  constitute  the 
state-of-the-art  in  LES-based  soot  prediction. 

3.1  Computational  Methodology  for  Chemical  Modeling  of  Soot 
Oxidation 

For  the  present  purpose  we  relied  on  theoretical  approaches  that  are  applicable  to  the  study 
of  large  systems. 

3.1.1  Energetics  of  Oxy radicals 

To  ensure  consistency  of  results,  geometry  optimization  of  oxyradicals  was  performed  at  two 
levels  of  theory.  One  was  the  generalized  gradient  approximation  (GGA)  with  the  Perdew- 
Burke-Ernzerhof  [33]  exchange-correlation  functional  and  double-zeta  quality  pseudo-atomic 
orbital  basis  set  [34]  as  implemented  in  the  SIESTA  software  [35].  The  force  tolerance 
in  conjugate  gradient  optimization  was  set  to  0.005  eV/AA.  The  other  approach  was  the 
PM6  semi-empirical  level  of  theory  [36]  in  its  restricted  open-shell  form  as  implemented  in 
MOPAC2009  software. 

The  harmonic  oscillator  model  of  aromaticity  (HOMA)  [37]  was  used  to  assess  local  elec¬ 
tronic  structure  of  individual  six-membered  rings  in  oxyradicals  in  order  to  explain  observed 
trends  in  relative  energies.  HOMA  is  a  geometry-based  criterion  that  connects  energetic  de¬ 
scriptions  of  aromaticity  to  the  geometry  of  the  system  that  does  not  become  prohibitively 
expensive  with  system  size.  The  HOMA  model  focuses  on  the  carbon-carbon  bond  lengths 
of  the  rings  contained  in  the  molecule,  and  is  defined  by, 

N 

HOMA=l--V(rcc-rcc  —benzene)  (l) 

n  z J 

2=1 

where  a  was  chosen  to  be  98.89  so  that  HOMA  =  0  for  a  Kekule  form  of  benzene  and 
HOMA  =  1  for  the  aromatic  form  of  benzene,  n  is  the  number  of  bond  lengths  in  the  ring, 
rcc  is  the  carbon-carbon  bond  length  in  the  system  under  consideration,  and  rcc-benzene  is 
set  to  the  experimental  carbon-carbon  bond  length  of  benzene.  Deviation  of  HOMA  from  1 
is  a  signature  of  deviation  of  the  aromatic  character  of  the  ring  from  that  of  benzene.  The 
sum  of  the  HOMA  values  for  individual  rings  in  the  PAH  (cumulative  HOMA)  can  be  used 
to  characterize  aromaticity  of  a  set  of  conjugated  rings.  In  PM6  calculations,  C-C  bonds  can 
be  shorter  than  the  values  used  to  parametrize  Eq.  1,  so  that  negative  values  of  HOMA  are 
encountered.  Technically,  this  indicates  a  ring  being  “less  aromatic”  than  a  Kekule  form  of 
benzene.  For  the  practical  purposes,  negative  values  of  HOMA  at  PM6  level  were  treated  as 
zero. 

It  should  be  noted  that  a  certain  ambiguity  is  encountered  in  optimization  of  oxyradical 
geometries.  Relaxation  of  the  carbon-carbon  framework  upon  oxidation  in  a  given  position 
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at  the  edge  may  result  in  several  alternative  arrangements  of  carbon-carbon  bonds  and, 
therefore,  different  isomers  with  respect  to  location  of  7r-aromatic  fragments.  It  is  impossible 
to  guarantee  that  the  most  stable  isomers  of  oxyradicals  are  found. 

3.1.2  Potential  Energy  Surfaces 

Density  functional  theory  (DFT)  was  employed  to  calculate  potential  energy  surfaces  of 
all  stable  species  and  transition  states  for  the  oxyradical  systems.  Geometry  optimization 
and  vibrational  frequency  calculations  were  performed  to  identify  all  stationary  points  on 
the  reaction  pathways  using  the  B3LYP  hybrid  functional  [38]  and  a  6-311G(d,p)  basis  set. 
Zero-point  energies  (ZPE)  and  vibrational  frequencies  were  scaled  by  a  factor  of  0.967  [39]. 
Transition  states  were  confirmed  to  connect  the  reactant  and  product  species  by  inspection  of 
the  normal  mode  for  the  single  imaginary  frequency  of  each  saddle  point.  All  the  quantum- 
chemical  calculations  were  carried  out  using  the  Gaussian  03  and  the  Gaussian  09  program 
packages. 

3.1.3  Reaction  Rate  Coefficients 

The  rate  coefficients  of  the  thermally  activated  reaction  systems  were  computed  using  the 
MultiWell  suite  of  codes  (MultiWell-2011.3)  [40].  MultiWell  solves  the  one-dimensional  time- 
dependent  energy-transfer  master  equations  for  a  multi-well  and  multi-channel  unimolecular 
reaction  system  using  the  Monte  Carlo  stochastic  method.  Microcanonical  rate  coefficients 
for  the  elementary  reactions  of  these  systems  were  calculated  with  MultiWell  at  the  Rice- 
Ramsperger-Kassel-Marcus  (RRKM)  level  of  theory.  This  approach  is  estimated  to  produce 
rate  coefficients  to  within  an  order-of-magnitude  accuracy. 

The  key  inputs  to  MultiWell — reaction  barriers,  frequencies,  and  moments  of  inertia — 
were  determined  from  the  DFT  calculations.  Following  Gilbert  and  Smith  [41],  the  real 
frequencies  below  150  cm-1  were  examined  by  graphically  visualizing  the  associated  normal 
mode  vibrations  to  identify  internal  rotational  modes.  All  internal  rotors  were  treated  as  1- 
D  hindered  rotors  with  rotational  potential  energy  barriers  calculated  at  the  B3LYP/3-21G 
level. 

Reaction  rates  were  computed  at  temperatures  ranging  from  1500  to  2500  K  and  pressures 
from  0.01  to  10  atm.  Argon  was  chosen  as  the  bath  gas  collider.  The  exponential-down  model 
with  (AEdown)  =  260  cm-1  was  used  to  describe  the  collisional  energy  transfer.  Lennard- 
Jones  parameters  were  estimated  from  an  empirical  correlation  [42],  The  exact  count,  with 
an  energy  grain  size  of  10  cm-1  for  the  first  segment  of  the  double  array  and  a  maximum 
energy  of  500,000  cm-1,  was  employed  to  determine  the  density  and  sum  of  states.  The 
maximum  energy  was  increased  to  2,000,000  cm-1  for  the  largest  11-ring  structure.  For  each 
set  of  initial  conditions,  the  number  of  stochastic  trials  was  varied  from  4  x  104  to  1  x  109 
to  keep  statistical  error  below  5  %. 

The  thermal  decomposition  rate  coefficients  were  derived  from  the  exponential  decay  of 
the  reactant  molecule  after  a  period  of  initial  relaxation,  as  outlined  and  tested  recently 
by  Golden  [43].  Following  this  procedure,  a  MultiWell  simulation  was  started  with  either 
a  monoenergetic  or  thermally  activated  initial  energy  distribution.  After  a  period  of  initial 
relaxation,  the  average  value  of  the  internal  energy  approached  a  constant  value,  and  a 


fraction  of  the  reactant  molecule  began  to  decompose  exponentially.  The  slope  of  the  decay 
on  a  semilog  plot  yielded  the  rate  coefficient  of  interest. 


3.2  LES  Model  for  Soot  Formation 

The  LES  soot  model  has  four  major  components:  a  statistical  model  for  evolution  of  the  soot 
particle  size  distribution  (Sec.  3.2.1),  a  gas-phase  combustion  model  (Sec.  3.2.2),  a  closure 
for  subfilter  correlations  of  soot  and  gas-phase  chemistry  (Sec.  3.2.3),  and  a  model  for  PAH 
mass  fraction  (Sec.  3.2.4).  Additionally,  Sec.  3.3  discusses  the  chemical  mechanism  used  for 
the  LES  simulations  presented  in  Sec.  4.3. 


3.2.1  Statistical  Soot  Modeling  with  HMOM 

Soot  particles  are  too  numerous  to  track  directly,  so  the  population  is  described  statistically 
with  the  Number  Density  Function  (NDF)  Nt.  Typically,  due  to  persistent  nucleation  of  soot 
particles  from  PAH,  the  NDF  is  bimodal  with  one  mode  containing  the  smaller  incipient 
spherical  particles  and  the  other  mode  containing  larger,  more  mature  fractal  aggregates 
composed  of  smaller  spherical  primary  particles.  In  order  to  describe  the  morphology  of 
these  aggregates,  the  NDF  is  defined  over  two  internal  dimensions:  the  volume  V  and  the 
surface  area  S.  Due  to  the  high-dimensionality  of  the  NDF  (three  spatial  dimensions  and 
two  internal  coordinates),  it  cannot  be  solved  for  directly,  and  a  statistical  model  is  required. 
For  application  to  LES,  due  to  cost  constraints,  the  only  tractable  model  is  the  Method  of 
Moments,  in  which  moments  of  the  NDF  are  solved  for  rather  than  the  NDF  itself.  The 
moments  of  the  NDF  of  order  x  and  y  in  volume  and  surface  area,  respectively,  Mxy  are 
given  by 

M.,,  =  F  V’S’N-  •  (2) 

i 

where  summation  over  i  implies  summation  over  all  particle  sizes  (written  here  in  a  dis¬ 
crete  sense).  While  the  Method  of  Moments  does  not  provide  the  distribution  directly,  the 
moments  of  the  NDF,  that  is,  the  soot  volume  fraction,  total  number  density,  average  pri¬ 
mary  particle  diameter,  etc.,  are  usually  the  quantities  of  interest,  particularly  in  engineering 
applications. 

The  evolution  of  the  NDF  is  governed  by  the  Population  Balance  Equation  (PBE),  and 
the  transport  equation  governing  the  evolution  of  the  moments,  obtained  by  taking  the 
moment  of  the  PBE,  is  given  by 


dMXjV  dujM^y 
dt  dxi 


d_ 

dxi 


v  dT 
Tdxi 


M, 


x,y 


M 


x,y  > 


(3) 


where  the  first  term  on  the  right  hand  side  is  the  thermophoresis  of  soot  particles.  Subse¬ 
quently,  this  term  will  be  combined  with  the  convective  term,  and  the  total  velocity  will  be 
denoted  by 


*  n  „vdT 

U'=U'~°-55fd^ 


(4) 


The  source  term  in  Eq.  3  contains  contributions  from  the  formation,  growth,  and  destruction 
processes  that  govern  the  internal  evolution  of  the  soot  population.  These  processes  include 
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particle  nucleation  from  PAH  dimers,  PAH  dimer  condensation,  particle  coagulation,  surface 
growth  by  the  hydrogen  abstraction/carbon  addition  (HACA)  mechanism,  oxidation,  and 
oxidation-induced  fragmentation.  Details  regarding  the  modeling  of  these  processes  can  be 
found  in  Mueller  et  al.  [18,  23,  25]  and  the  references  therein. 

The  major  challenge  with  Method  of  Moments  is  that,  in  general,  the  source  terms  are 
unclosed.  In  other  words,  the  source  terms  for  the  moment  equations  depend  on  moments 
which  are  not  solved  for.  In  this  work,  closure  is  obtained  with  the  Hybrid  Method  of 
Moments  (HMOM)  [23].  Briefly,  in  HMOM,  the  contribution  to  the  moments  from  the 
smaller  incipient  particles  is  described  with  a  delta  function  as  in  the  Direct  Quadrature 
Method  of  Moments  (DQMOM)  [44] ,  and  the  contribution  to  the  moments  from  the  larger 
particles  is  described  with  polynomial  interpolation  as  in  the  Method  of  Moments  with 
Interpolative  Closure  (MOMIC)  [20].  An  arbitrary  moment  is  then  given  by 


M, 


x,y 


=  V?SJJVo 


R  r 

exP  ^2  ^2  ar,kxkyr~k 

r= 0  k= 0 


(5) 


where  the  location  of  the  delta  function  (Vo,  5o)  is  fixed  and  R  is  the  order  of  the  polynomial 
interpolation.  Given  the  weight  of  the  delta  function  Nq,  the  polynomial  interpolation  coeffi¬ 
cients  are  obtained  from  a  set  of  moments  which  are  solved  for.  In  this  work,  R  —  1,  so  three 
moment  equations  are  solved:  the  total  number  density  Afore  the  total  soot  volume  Mqo,  and 
the  total  soot  surface  area  Af0,i-  In  diffusion  flames,  the  primary  soot  growth  mechanisms 
are  nucleation  and  condensation  [24],  and  the  soot  volume  fraction  is  not  sensitive  to  the 
number  of  moments.  In  addition,  a  transport  equation  is  solved  for  the  weight  of  the  delta 
function  A0.  The  transport  equation  for  this  quantity  is  the  same  as  Eq.  3  with  the  source 
term  given  by 


M-a- 


0 


No  =  lim  „ 

a,0—>co  v0~asrp 


(6) 


'0 


Additional  details  regarding  HMOM  including  the  modeling  of  the  formation,  growth,  and 
destruction  processes  can  be  found  in  Mueller  et  al  [23,  25]. 


3.2.2  Combustion  Model 

For  sooting  flames,  the  combustion  model  must  be  able  to  predict  accurately  not  only  the 
oxidation  of  the  fuel  but  also  the  formation  of  soot  precursors  such  as  acetylene  and  PAH. 
However,  at  the  same  time,  the  model  must  be  tractable  in  terms  of  both  cost  and  closure. 
Therefore,  the  thermochemical  state  will  be  described  with  the  Flamelet/Progress  Variable 
(FPV)  model  of  Pierce  and  Moin  [45].  In  the  FPV  model,  the  local  thermochemical  equa¬ 
tion  of  state  is  obtained  from  the  solutions  of  the  steady  flamelet  equations  [46],  which  are 
parameterized  by  the  mixture  fraction  Z  and  a  reaction  progress  variable  C : 

t  =  F(Z,C),  (7) 

where  £  is  any  thermochemical  variable  (density,  temperature,  mass  fractions,  etc.)  and  T 
is  the  functional  relationship  obtained  from  the  solution  of  the  steady  flamelet  equations. 

Due  to  the  long  time  scales  on  which  soot  evolves  combined  with  the  sensitivity  of  soot 
evolution  to  temperature  (and  the  thermochemical  state  in  general),  radiation  losses  must  be 
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accounted  for.  Ihme  and  Pitsch  [47]  extended  the  FPV  model  to  account  for  heat  losses  due  to 
radiation  (RFPV).  The  steady  flamelet  database  is  augmented  with  specific  solutions  of  the 
unsteady  flamelet  equations,  and  a  heat  loss  parameter  H  is  added  to  the  parameterization. 
The  thermochemical  equation  of  state  is  then  given  by 

Z  =  g(Z,C,H),  (8) 

where  Q  is  the  functional  relationship  obtained  from  the  augmented  flamelet  database. 

For  sooting  flames,  some  quantities,  such  as  the  source  terms  in  Eq.  3,  are  functions  of 
the  both  thermochemical  quantities  and  the  soot  scalars.  Therefore,  the  equation  of  state  is 
generalized  to  allow  for  a  dependence  on  both  sets  of  quantities: 

<t>  =  J{G(Z,C,H),Mi)  ,  (9) 

where  (f)  is  now  any  quantity  that  could  depend  on  not  only  the  thermochemical  state  but 
also  the  soot  scalars  and  AF,  is  the  vector  of  soot  scalars,  that  is,  the  moments  Mx  y  and 
the  weight  of  the  delta  function  N0.  For  all  of  the  quantities  needed  in  the  soot  model,  the 
general  equation  of  state  J  can  be  written  as  the  product  of  a  function  that  depends  only 
on  the  thermochemical  state  and  one  that  depends  only  on  the  soot  scalars: 

cf)  =  g{Z,C,H)lC{Mi)  ,  (10) 


where  K,  is  unity  if  0  is  a  strictly  thermochemical  quantity. 

A  convenient  definition  of  the  mixture  fraction  is  a  conserved  scalar  that  is  defined  to 
be  zero  in  the  oxidizer  stream  and  one  in  the  fuel  stream  [48].  A  major  benefit  of  this 
definition  is  that  the  mixture  fraction  transport  equation  is  independent  of  the  transport 
model  used  for  individual  species.  However,  in  sooting  flames,  the  conserved  scalar  definition 
is  inappropriate.  In  order  to  form  soot,  PAH  are  removed  from  the  gas-phase,  and  the  mixture 
is  locally  leaned.  Therefore,  the  mixture  fraction  will  be  defined  by  the  following  transport 
equation: 


dpZ 

dt 


dpuiZ 

dxi 


d 

dxi 


+  rnz  , 


(11) 


where  the  mixture  fraction  Lewis  number  Le^  is  unity.  The  mixture  fraction  source  term 
rhz  is  defined  analogously  to  the  element  mass  fraction  based  mixture  fraction  of  Bilger  [49] : 
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(12) 


where  m,  are  element  mass  source  terms,  p  is  the  source  term  in  the  continuity  equation 
due  to  the  removal  of  PAH  from  the  gas-phase,  Z,  are  the  element  mass  fractions  with  a 
second  subscript  indicating  the  value  in  the  fuel  or  oxidizer  stream,  zy  are  the  stoichiometric 
coefficients,  and  H  j  are  the  molar  masses.  With  the  source  term  defined  in  this  way,  in 
the  limit  of  unity  Lewis  numbers  for  all  species,  this  mixture  fraction  definition  would  be 
equivalent  to  the  element  mass  fraction  based  mixture  fraction  of  Bilger  [49].  With  the 
mixture  fraction  defined  as  in  Eq.  11,  additional  terms  appear  in  the  flamelet  equations. 
These  will  be  discussed  below. 
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In  previous  works  utilizing  the  FPV  model  [45,  47],  the  progress  variable  is  defined  as 
the  sum  of  the  mass  fractions  of  the  major  products  of  combustion.  However,  in  sooting 
flames,  this  conventional  definition  is  again  not  appropriate.  The  PAH  that  are  removed 
from  the  gas- phase  have  a  much  smaller  H/C  ratio  than  typical  aliphatic  fuels  (less  than  one 
compared  to  greater  than  two).  As  a  result  of  removing  PAH  from  the  gas-phase,  the  local 
H/C  ratio  increases,  and  the  local  effective  fuel  composition  is  changed.  To  account  for  this 
effect,  the  progress  variable  is  defined  by  the  following  transport  equation: 

dpC  dpUiC  _  d  f  dC\  mSYc 
dt  +  dxi  ~  dxi  \f  CdxJ  +  C* 

where  the  progress  variable  Lewis  number  Lee  is  unity,  msyc  is  the  conventional  progress 
variable  source  term,  that  is,  the  sum  of  the  source  terms  of  the  major  productions  of 
combustion,  and  C*  is  a  normalizing  factor  that  accounts  for  the  change  in  stoichiometry 
due  to  the  removal  of  PAH.  In  addition,  the  progress  variable  is  defined  to  be  zero  in  both  the 
fuel  and  oxidizer  streams,  irrespective  of  composition.  The  normalizing  factor  C*  is  defined 
as  the  sum  of  the  mass  fractions  of  CO2  and  H20  based  on  an  idealized  one-step  reaction 


CaHft  +  c02  +  dN2  — »  eC02  +  / H20  +  c/02  +  hN2  ,  (14) 


where  the  fuel  composition  and  stoichiometric  coefficients  are  determined  from  the  local  H/C 
ratio,  a  constant  reference  (element  mass  fraction  based)  mixture  fraction  at  the  nominal 
H/C  ratio,  and  the  local  N/O  ratio.  As  the  local  H/C  ratio  increases,  c  increases,  and, 
therefore,  C*  decreases.  The  progress  variable  defined  in  this  way  ensures  that  the  equation 
of  state  (Eqs.  7-10)  is  unique.  In  non-sooting  flames,  the  H/C  ratio  is  constant;  C*  is 
therefore  a  constant;  and  the  progress  variable  as  defined  here  is  a  linear  function  of  the 
conventional  progress  variable. 

Finally,  the  role  of  the  third  coordinate  in  the  flamelet  tabulation  is  to  distinguish  between 
states  on  which  radiation  has  had  little  effect  and  states  for  which  radiation  has  lowered  the 
enthalpy.  One  possible  choice  for  this  coordinate  is  the  enthalpy  itself  or  the  corresponding 
enthalpy  deficit.  The  latter  quantity  has  the  benefit  of  a  well  defined  adiabatic  boundary 
(zero)  for  the  simplest  cases.  However,  when  a  detailed  transport  model  is  used  or  when 
soot  is  present  (affects  the  adiabatic  boundary  for  the  enthalpy  deficit),  a  transport  equation 
definition  is  more  convenient.  This  quantity,  which  will  be  called  the  heat  loss  parameter  H , 
is  defined  by  the  transport  equation 


dpH  dpu,H 
dt  dxi 


IPO 


+  pH  +  ()rad  > 


(15) 


where  the  heat  loss  parameter  Lewis  number  Le#  is  unity  and  ({rad  is  the  radiation  source 
term.  The  heat  loss  parameter  H  is  defined  to  be  initially  zero  and  becomes  increasingly 
negative  as  the  cumulative  effect  of  radiation  becomes  stronger.  The  second  term  on  the 
right  hand  side  ensures  that  H  is  identically  zero  everywhere  if  the  radiation  source  term 
is  zero.  The  radiation  model  used  in  this  work  is  an  optically  thin  gray  gas  approach  [50]. 
In  the  validation  case  considered  in  the  next  section,  the  soot  volume  fraction  is  sufficiently 
small  that  soot  radiation  can  be  neglected. 
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As  mentioned  previously,  the  flamelet  equations  contain  additional  terms  other  than  those 
found  in  Peters  [46]  due  to  the  source  terms  in  the  mixture  fraction  and  continuity  equations 
For  unity  Lewis  numbers,  which  is  the  transport  model  considered  in  this  work,  the  flamelet 
equations  for  the  species  mass  fractions  are 

The  final  two  terms  on  the  right  hand  side  are  the  new  terms  consisting  of  a  source  term  due 
to  the  density  source  term  and  a  convective  term  in  mixture  fraction  space  due  to  both  the 
density  and  mixture  fraction  source  terms.  The  primary  role  of  these  terms  is  to  ensure  that 
mass  is  conserved,  and  they  do  not  play  a  substantial  role  in  the  dynamics  of  the  flamelet 
equations.  The  flamelet  equation  for  the  temperature  is 


dT  pcpXd2T  hPXdTdcp+^2 


pCp  dr  ~ 


2  dZ2  2  dZdZ 


pCpiX  dYi  dT  ■  ~  .  dT  .  dT 

b-H-2^ himi+pcPZ—-cpmz  —  +qR AD  , 


2  dZdZ 


(17) 

where  H  is  the  enthalpy  source  term  due  to  the  removal  of  PAH  from  the  gas-phase.  In 
addition  to  these  quantities,  flamelet  equations  are  also  solved  for  the  progress  variable  C 
and  the  heat  loss  parameter  H.  From  the  transport  equation  for  C  (Eq.  13),  the  flamelet 
equation  is 

dC  px  d2C  rhxyc  .  ,  8C 


9  Bt  2  dZ 2 


C* 


~  pC  +  (pZ  -  mz)  , 


and,  from  the  transport  equation  for  H  (Eq.  15),  the  flamelet  equation  is 

dH  pX  d2H  dH 

P—  =  — VV7  +  (Pz  -mz)-7^  +  <?rad  • 


dr 


2  dZ2 


(18) 


(19) 


The  values  of  C  and  H  obtained  from  the  solutions  of  these  flamelet  equations  are  then  used 
to  parametrize  the  solutions  of  Eqs.  16,  17  with  the  RFPV  combustion  model. 


3.2.3  Closure  of  the  joint  subfilter  PDF 

In  LES,  spatially  Favre  filtered  quantities  are  evolved  for  the  velocity  and  scalars  (Z,  C , 
H,  and  At,),  and  spatially  Favre  filtered  quantities  are  needed  from  the  equation  of  state 
to  evolve  these  equations  (density,  source  terms,  etc.).  These  Favre  filtered  quantities  are 
obtained  from  the  equation  of  state  by  convolution  with  the  joint  subfilter  PDF  P: 


<KZ,  C,H,Mi ) 


J  (Z,  C,  H,  Mi)  P  (Z,  C,  H,  Mi)  dMidHdCdZ  .  (20) 


Two  basic  options  exist  for  specifying  this  joint  PDF,  the  presumed  shape  PDF  approach 
and  the  transported  PDF  approach,  both  of  which  are  considered  in  this  project. 


3.2.3. 1  Presumed  Shape  PDF  Approach  The  model  for  the  joint  subfilter  PDF  is 
developed  and  validated  a  priori  by  Mueller  and  Pitsch  [30]  for  any  general  combustion 
model.  In  the  remainder  of  this  section,  a  brief  overview  of  this  model  is  given  including 
details  specific  to  the  RFPV  combustion  model  considered. 
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As  a  first  step  in  modeling  the  joint  subfilter  PDF,  Eq.  10  and  conditional  distributions 
are  used  to  split  the  subfilter  PDF  into  a  thermochemical  component  and  a  soot  distribution 
conditioned  on  the  thermochemical  variables: 


(f>(Z,  C,H,Mi ) 


g  (Z,  C,  H )  K  (Mi)  P  (Z,  C,  H )  P  ( Mi\Z ,  C,  H )  dMidHdCdZ  . 

(21) 


Furthermore,  due  to  the  disparate  time  scales  between  the  evolution  of  the  thermochemical 
state  and  soot,  Mueller  and  Pitsch  [30]  argue  that  the  conditional  soot  distribution  can  be 
modeled  by  the  marginal  distribution.  The  previous  relationship  then  simplifies  to 


<KZ,C,  H,  Mi) 


g  (Z,  C,  H)  P  (Z,  C,  H)  dHdCdZ 


K(Mi)P(Ml)dMi ,  (22) 


and  the  thermochemical  and  soot  portions  of  the  convolution  are  now  completely  indepen¬ 
dent.  The  model  for  the  joint  subfilter  PDF  in  Eq.  20  has  now  been  reduced  to  modeling 
the  thermochemical  subfilter  PDF  and  the  marginal  soot  subfilter  PDF. 

The  thermochemical  subfilter  PDF  for  the  RFPV  model  is  from  the  work  of  Ihme  and 
and  Pitsch  [47].  Mathematically,  the  thermochemical  PDF  is  expressed  by  first  recasting 
the  equation  of  state  (Eq.  8)  in  terms  of  two  quantities  that  uniquely  identify  each  flamelet 
solution  in  the  database:  A  =  C  ( Zst )  and  $  =  H  (Zst).  Thermochemical  quantities  are  then 
given  by 

Z  =  g(z,c,H)  =  g*(z,  a,$)  (23) 

or,  as  a  filtered  quantity, 


£  = 


g  (Z,  C,  H)  P  (Z,  C,  H)  dHdCdZ  = 


g *  (Z,  A,  <F)  P  (Z,  A,  <F)  d<$>d\dZ  . 

(24) 

By  definition,  Z,  A,  and  $  are  independent,  and  the  marginal  distributions  are  modeled  as 
a  beta  distribution  for  the  mixture  fraction  [51-53]  and  delta  distributions  for  A  and  <F  [47]: 

P  (Z,  A,  d>)  =  fZ;  Z,  Z S  (A  -  X)  .  (25) 


The  resulting  dependence  on  A  and  <F  is  then  reexpressed  as  a  dependence  on  C  and  H, 
assuming  a  unique  inversion  exists,  which  is  guaranteed  if  Eq.  8  is  unique.  In  words,  this 
mathematical  formalism  is  practically  implemented  by  individually  convoluting  each  flamelet 
solution  in  the  database  with  the  beta  distribution  for  the  mixture  fraction  and  tabulating 
the  resulting  set  of  data  as  a  function  of  the  filtered  mixture  fraction,  subfilter  mixture 
fraction  variance,  filtered  progress  variable,  and  filtered  heat  loss  parameter. 

The  transport  equation  for  the  filtered  mixture  fraction  is  obtained  by  spatially  filtering 
Eq.  11  and  is  given  by 


dpZ  dpUiZ 
dt  dxi 


_d_ 

dxi 


pUiZ  —  puiZ 


d_ 

dxi 


+  rhz  ■ 


(26) 


The  filtered  source  term  has  been  closed  with  the  use  of  the  presumed  subfilter  PDF  (Eq.  25). 
Likewise,  the  transport  equation  for  the  filtered  progress  variable  is  given  by 


dpC  dpUiC  d 

dt  dxi  dxi 


ripYf  \ 
C*  )  ’ 


(27) 
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and  the  transport  equation  for  the  filtered  heat  loss  parameter  is  given  by 


dpH  dptiiH 
dt  dxi 


A 

dxi 


pUiH  —  pUiH  )  + 


d_ 

dxi 


+  pH  +  gRAD  . 


(28) 


For  the  mixture  fraction  used  in  this  work,  an  algebraic  model  for  the  corresponding  sub¬ 
filter  variance  could  be  derived  by  writing  the  transport  equation  for  the  subfilter  variance, 
assuming  local  equilibrium,  and  replacing  the  subfilter  scalar  dissipation  rate  with  a  linear 
relaxation  model  [54],  The  modeled  subfilter  variance  would  then  be  given  by 


Z"2 


Cy 

pS 


dZ  dZ 


dxi  dxi 


( pZ 2  -  pZ 2) 


+  2  mzZ  —  2mzZ 


(29) 


where  S  is  the  magnitude  of  the  filtered  strain  rate  tensor,  Dt  is  the  subfilter  diffusivity 
used  to  close  the  scalar  flux  in  the  mixture  fraction  equation,  and  Cy  is  a  constant  that 
must  specified  or  determined  dynamically.  The  last  two  terms  on  the  right  hand  side  are 
due  to  the  source  terms  that  appear  in  the  continuity  and  mixture  fraction  equations  and 
are  themselves  a  function  of  the  subfilter  variance.  As  a  result  of  these  terms,  a  dynamic 
procedure  for  Cy  would  require  an  iterative  approach,  with  test  filtering  required  at  each 
iteration.  Depending  on  the  number  of  iterations  needed  for  convergence,  this  model  could 
become  very  expensive  and  is  therefore  discarded  as  computationally  impractical.  Instead, 
a  simpler  and  less  expensive  approach  is  taken,  and  the  subfilter  mixture  fraction  variance 
is  obtained  from  the  solution  of  a  transport  equation  for  the  filtered  square  of  the  mixture 
fraction  ( Z "2  =  Z2  —  Z2): 


dpZ 2  dpUiZ2 
dt  dxi 


—2pDz 


dZ  dZ 

dxi  dxi 


-pXsgs-pZ2+2mzZ  , 

(30) 


where  ysgs  is  the  subfilter  contribution  to  the  dissipation  rate.  This  equation  is  preferred 
over  a  transport  equation  for  the  subfilter  variance  directly  since  a  production  term,  which 
is  susceptible  to  large  numerical  errors,  is  not  computed  explicitly  [55,  56].  Following  the 
work  of  Ihme  and  Pitsch  [57] ,  the  subfilter  scalar  dissipation  rate  is  modeled  with  a  linear 
relaxation  model  with  the  timescale  determined  from  the  subfilter  eddy  viscosity: 


w  =  cx^  (z2  -  Z2) 


(31) 


where  Cx  =  20. 

Due  to  the  large  spatial  intermittency  of  soot,  Mueller  and  Pitsch  [30]  modeled  the  soot 
subfilter  PDF  as  a  double  delta  distribution  with  a  ” non-sooting”  mode  and  a  ’’sooting” 
mode: 

P  (Mi)  =  cvS  (Mi)  d-(l  -u)8  (Mi  -  MX)  ,  (32) 

where  u>  is  the  subfilter  intermittency  and  the  M*  are  defined  such  that  convolution  with 
the  PDF  recovers  the  filtered  values  of  the  soot  moments  Mi.  The  subfilter  intermittency  is 
smaller  when  soot  is  more  homogeneous  at  subfilter  scales  and  can  be  determined  from  the 
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filtered  square  of  any  one  of  the  soot  moments.  However,  Mueller  and  Pitsch  [30]  showed 
that  the  model  performs  best  if  the  number  density  is  used,  and  the  subfilter  intermittency 
is  given  by 

M0,o2 

UJ  =  1  —  =!=  . 

M0,o2 

The  filtered  square  of  the  number  density  is  obtained  from  the  solution  of  a  transport  equa¬ 
tion. 

The  transport  equation  for  the  filtered  soot  moments  (and  weight  of  the  delta  function 
N0)  is  obtained  by  applying  a  spatial  filter  to  Eq.  3  and  is  given  by 


dMx^y  dtijMxg 
dt  dxi 


_d_ 

dxi 


+  M 


x,y  i 


(34) 


and  the  corresponding  equation  for  the  filtered  square  of  the  number  density  is  given  by 


dMo,o2  chq*Mo,o2 
dt  dxi 


A 

dxi 


w^du- 

OXi 


+  2Mo;oMo,o  • 


(35) 


In  the  latter  equation,  the  subfilter  contribution  to  the  product  of  the  number  density  squared 
and  the  divergence  of  the  total  velocity  has  been  neglected  following  Mueller  and  Pitsch  [30] . 
Due  to  the  simple  nature  of  the  soot  subfilter  PDF,  convolution  with  the  double  delta 
distribution  occurs  on-the-fly,  and  the  portions  of  the  soot  source  terms  which  are  functions 
of  the  thermochemical  state  are  stored  in  the  same  table  as  for  the  other  quantities  needed 
for  the  RFPV  combustion  model. 


3. 2. 3. 2  Transported  PDF  approach  In  the  transported  PDF  approach,  a  partial  dif¬ 
ferential  equation  describing  the  evolution  of  the  PDF  in  physical  coordinates  and  the  ther¬ 
mochemical  compositional  space  —  (X-  6',  H,  At,)  is  directly  solved.  In  this  study,  the 
transported  PDF  approach  is  employed  for  simulation  of  a  turbulent  jet  flame  in  cylindrical 
coordinates.  For  this  purpose,  the  transport  equation  is  written  in  cylindrical  coordinates 
as  follows  [58]: 


OGl  d_ 
dt  dr 


Ar  -\ - )  Gl 

r 


vt\4>  +  S„(0) )  ^ 


+  T  (ACGL)  +  T  (BGL)  +  T  (®Gi 

d2  d 

+  d^(yBGL">  =  ~d^a 

where  G/  =  rP.  B  =  Dt-  Ar,Ag,  and  Az  are  given  by 

1  d 

Ar  =  ur  +  ( pB ) 

^  =  +  Tplie (m 

a  ~  1  d 

A‘=u>+ wDB)- 


(36) 


(37) 

(38) 

(39) 
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The  terms  on  the  RHS  of  the  PDF  transport  equation  represent  transport  in  composition 
space  and  consist  of  molecular  mixing  and  chemical  reactions,  respectively.  The  primary 
advantage  of  the  PDF  approach  is  that  the  chemical  source  term  appears  closed  and  does 
not  require  modeling.  However,  modeling  is  required  for  describing  the  conditional  diffusion 

term,  A4i\(p: 

T>i\4>  =\\7  ■  pDS74>i\4>.  (40) 

P 

This  term  is  closed  using  the  interaction  by  exchange  with  the  mean  (IEM)  model  [59-61]: 

A  (aw)  =  v  •  ?5va/P + ±  [i  (a  -  a)  a]  ■  («) 


Here  r  is  a  mixing  time- scale  that  needs  to  be  specified,  and  D  is  a  common  diffusivity  that 
is  used  for  all  the  scalars  in  ip.  In  this  work,  this  is  taken  to  be  equal  to  the  diffusivity 
of  mixture  fraction.  Note  that  the  differential  diffusion  between  soot  moments  and  the 
gas  phase  scalars  could  be  important  [62],  As  a  first  step,  this  is  neglected  here,  but  the 
differential  diffusion  of  high-Schmidt  number  scalars  could  be  considered  by  modifying  the 
dissipation  time  scale  [63,  64], 

The  transported  PDF  approach  requires  solution  of  the  high- dimensional  PDF  transport 
equation  (Eq.  36),  which  spans  eleven  dimensions  for  the  state  vector  used  in  this  work. 
Conventional  finite  difference  or  finite  volume  methods  are  not  tractable,  and  a  Lagrangian 
Monte-Carlo  method  is  used  [59-61,  65,  66].  In  the  Lagrangian  approach,  the  computational 
domain  is  seeded  with  a  large  number  of  notional  particles  that  evolve  in  physical  and 
compositional  spaces  using  a  set  of  stochastic  differential  equations.  Each  of  these  particles 
carries  a  property  vector  that  consists  of  the  particle  location  vector,  a  characteristic  weight, 
and  the  ip  vector. 

The  evolution  of  the  particles  in  physical  space  is  governed  by  the  following  equations: 


drn 

d9n 

dzn 


(42) 

(43) 

(44) 


where  the  superscript  n  denotes  the  particle  index,  and  dW%  is  a  Weiner  diffusion  process 
with  zero  mean  and  variance  of  1,  and  B\  =  D  +  Dt .  The  velocity  and  diffusivity  fields 
used  in  this  equation  are  interpolated  onto  the  particle  locations  using  trilinear  interpolation. 
The  transport  in  composition  space  consists  of  mixing  and  chemical  reactions.  The  particle 
equation  corresponding  to  the  IEM  mixing  model  is  given  by 

d(pn  =  *  (y  ~  </>n)  ,  (45) 


where  ip  is  the  filtered  composition  vector  in  a  given  filter  volume,  <pn  is  the  particle  compo¬ 
sition  vector,  and  r  is  a  mixing  time  scale.  In  this  work,  the  mixing  time  scale  is  determined 
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using  the  local  effective  diffusivity  and  the  filter  scale. 


A2 


T  =  Cfi 


D  +  Dt ’ 


(46) 


where  c p  is  a  model  coefficient.  The  dynamic  procedure  proposed  by  Raman  and  Pitsch  [66] 
is  used  for  evaluating  c^.  The  particle  equations  thus  reduce  to  solving  a  vector  of  ordinary 
differential  equations. 

d<pn  =  S  (</>n)df,  (47) 

where  S  is  the  vector  of  chemistry  and  soot  source  terms. 

The  particle  weights  are  proportional  to  the  mass  of  the  fluid  they  represent.  It  should  be 
noted  that  the  PDF  transport  equation  is  formulated  for  rFL  instead  of  FL.  This  modified 
formulation  ensures  that  the  weights  have  no  evolution  equation.  In  other  words,  the  particle 
weights  are  initialized  at  the  inlet  but  do  not  change  with  time.  To  ensure  consistency,  the 
particle  weights  are  set  proportional  to  pnr^  at  the  inlet,  where  p"  is  the  Eulerian  density 
interpolated  to  the  particle  location,  and  rfi  is  the  radial  location  of  the  particle. 


3.2.4  PAH  Model 

The  final  component  of  the  LES  model  to  be  discussed  is  a  model  for  obtaining  the  PAH  mass 
fraction  and,  subsequently,  the  mass  transfer  rate  from  the  gas-phase  to  soot.  The  simplest 
approach  would  be  to  take  the  PAH  mass  fraction  directly  from  the  RFPV  combustion 
model.  However,  the  DNS  study  of  Bisetti  et  al.  (see  Ref.  24  and  Sec.  4.2)  demonstrated 
that  the  PAH  mass  fraction  deviates  significantly  from  the  steady  flamelet  model  due  to  the 
relatively  slow  PAH  chemistry.  As  the  scalar  dissipation  rate  changes,  species  with  faster 
chemistry  (e.g.,  the  major  species  of  combustion)  respond  quickly  to  these  changes,  and  the 
unsteady  term  in  the  flamelet  equations  is  negligible.  For  species  with  slower  chemistry,  a 
longer  time  is  needed  to  respond  to  changes  in  the  dissipation  rate;  the  unsteady  term  in 
the  flamelet  equations  must  be  retained;  and  the  mass  fractions  of  these  species  cannot  be 
accurately  described  with  the  solutions  to  the  steady  flamelet  equations. 

Like  PAH,  NO  is  governed  by  relatively  slow  chemistry.  To  model  the  unsteady  effects, 
Ihme  and  Pitsch  [47]  developed  a  transport  equation  model  for  NO.  The  same  approach  will 
be  taken  here  for  PAH.  The  spatially  filtered  transport  equation  for  a  lumped  PAH  mass 
fraction  is  given  by 


<9pEpAH  dpUiYp ah  0  .  __  _  , , 

-m-  +  ^ PAH  -  wipAH 


+  w, 1  pDfah 


dYpAH  \ 


dxi 


1 


+  fflPAH  J  (48) 


where  rhpAH  is  the  sum  of  the  sources  of  all  of  the  PAH  species.  If  all  of  the  PAH  have  the 
same  molecular  diffusivity  Dpah  (e.g.,  a  unity  Lewis  number  assumption),  then  the  preceding 
equation  is  exactly  the  transport  equation  for  the  sum  of  the  PAH  mass  fractions.  If  the  PAH 
do  not  all  have  the  same  molecular  diffusivity,  then  the  diffusivity  of  naphthalene  (generally 
the  PAH  species  with  the  largest  concentration)  is  preferred  for  DPah5  and  the  preceding 
equation  is  then  an  approximate  model  for  a  lumped  PAH. 

With  a  conventional  model  for  the  scalar  flux,  the  lone  unclosed  term  in  Eq.  48  is  the 
filtered  source  term.  This  source  can  be  decomposed  into  three  components:  a  chemical 
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production  term  rh+,  a  chemical  consumption  term  m_ ,  and  an  additional  consumption 
term  due  to  the  mass  transfer  rate  from  the  gas-phase  to  soot  (dimerization)  mo-  For  each 
individual  species,  the  chemical  production  term  is  independent  of  the  concentration  of  the 
species;  the  chemical  consumption  term  is  linear  in  the  concentration  of  the  species;  and  the 
dimerization  term  is  quadratic  in  the  concentration  of  the  species.  For  the  lumped  species, 
the  source  term  will  be  decomposed  in  the  same  way: 


hr  pah  =  (rh+)  + 


FpAH  + 


(49) 


Written  in  this  way,  the  three  terms  in  parentheses  are  nearly  independent  of  the  PAH  mass 
fraction  and  can  be  taken  directly  from  the  RFPV  combustion  model.  Following  the  work 
of  Ihme  and  Pitsch  for  NO  [47],  the  correlations  in  the  second  and  third  terms  are  closed  by 
assuming  scale-similarity  between  the  transport  equation  model  and  the  RFPV  combustion 
model.  The  final  form  of  the  source  term  is  then 


—  —RFPV  —RFPV 

mpAH  =  rn+  +  m_ 


—RFPV 

+  mD 


(  VpAH  A 


(50) 


where  the  superscript  RFPV  indicates  quantities  which  are  obtained  from  the  RFPV  combus¬ 
tion  model. 


3.3  Chemical  Mechanism 

For  construction  of  the  flamelet  look-up  tables  used  in  the  LES  computations,  solutions  to 
the  steady  and  unsteady  flamelet  equations  are  obtained  with  FlameMaster  [67].  The  base 
chemical  mechanism  is  from  Blanquart  et  al.  [68]  with  extensions  from  Narayanaswamy  et 
al.  [69].  The  mechanism  has  been  extensively  validated  for  simple  C0-C4  fuels,  n-heptane, 
iso-octane,  benzene,  and  several  substituted  aromatics  (toluene,  ethylbenzene,  m-xylene, 
and  «-methylnaphthalene)  at  both  lean  and  rich  conditions  in  several  configurations.  In 
addition,  the  mechanism  contains  PAH  chemistry  up  to  Ci6  and  Ci«  species  (pyrene,  etc.). 
Numerous  PAH  formation  and  growth  pathways  are  included  in  the  mechanism. 

For  the  commercial  natural  gas  mixture  considered  in  this  work,  the  mechanism  predicted 
that  most  naphthalene  was  formed  from  the  combination  of  an  indenyl  radical  and  a  methyl 
radical  [70,  71].  However,  the  as-compiled  mechanism  did  not  contain  a  pathway  for  the 
decomposition  of  indenyl  under  rich  conditions,  and  all  of  the  indenyl  that  was  formed  (pri¬ 
marily  through  acetylene  addition  to  a  benzyl  radical)  eventually  ended  up  as  naphthalene, 
leading  to  unphysically  large  PAH  concentrations.  Therefore,  a  global  rate  for  the  thermal 
decomposition  of  indenyl  was  added  to  the  mechanism  based  on  the  work  of  Laskin  and 
Lifshitz  [71]: 

C9H7  C3H3  +  C4H2  +  C2H2  (51a) 

C9H7  C5H5  +  C4H2  (51b) 

where  the  branching  ratio  between  the  two  product  channels  is  two  in  favor  of  the  first. 
At  lean  conditions,  this  thermal  decomposition  is  not  competitive  with  oxygen- containing 
species  and,  therefore,  does  not  affect  any  of  the  previous  validation  results  [68,  69]. 
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4  Results  and  Discussion 


This  section  begins  by  presenting  results  from  our  investigation  of  soot  oxidation  (Sec.  4.1), 
then  proceeds  to  discuss  a  series  of  turbulent  flame  investigations,  including  a  DNS  study 
of  soot  formation  (Sec.  4.2)  and  validation  studies  of  the  LES  soot  model,  using  both  the 
presumed  PDF  (Sec.  4.3.3)  and  transported  PDF  (Sec.  4.3.4)  approaches. 

4.1  Oxidation  of  Soot 

We  begin  this  section  by  first  reporting  results  of  our  study  of  bay-capping  reactions,  one 
of  the  important  reaction  steps  in  soot  surface  growth,  and  examine  their  effect  on  the 
evolution  of  soot  particle  growth.  We  then  examine  the  thermodynamic  stability  of  the 
key  intermediates  of  soot  oxidation,  graphene-edge  oxyradicals,  located  at  both  zigzag  and 
armchair  sites,  and  then  proceed  with  the  kinetics  of  graphene  oxyradical  decomposition. 

4.1.1  Surface  growth:  Bay-capping  reactions 

Controlling  formation  of  soot  in  modern  combustion  systems  requires  a  good  mechanistic 
understanding  of  soot  particle  growth.  One  of  the  biggest  unknowns,  particle  inception, 
received  substantial  attention  in  recent  years  and  its  mechanistic  picture  steadily  advances 
toward  a  unified  theory  [5] .  Another  part  of  the  mechanism,  surface  growth,  has  transformed 
from  an  empirical  rate  law  to  a  simple  elementary-reaction  model,  to  a  complex,  multi-step 
process.  The  complexity  of  the  surface  growth  mechanism  stems  from  numerous  possibilities 
of  chemical  transformations  taking  place  at  graphene  edges,  presumed  to  be  the  building 
units  of  soot  particles.  Among  them  are  game-changing  processes:  migrations  of  lone  [72] 
and  embedded  [73]  five-member  rings,  which  give  rise  to  rapid  lateral  translation  of  five- 
member  rings  along  zigzag  edges. 

Most  recently,  Whitesides  and  Frenklach  [10]  presented  a  new  detailed  kinetic  Monte- 
Carlo  (KMC)  model  of  graphene-edge  growth  with  a  total  of  42  surface  transformations. 
Evolving  edge  morphology  and  growth  rate  were  found  to  be  affected  by  the  rates  of  key 
reactions.  One  of  the  most  interesting  among  them  is  capping,  the  addition  of  acetylene  to 
an  embedded  five- member  ring,  shown  in  Fig.  1.  The  significance  of  this  reaction  step  is 


+  H 


Figure  1:  Bay-capping  reaction. 


several-fold  [10].  First,  it  pins  down  the  five-member  ring  and  prevents  it  from  migrating. 
As  ring  migration  generally  leads  to  smoother  surfaces,  this  aspect  of  capping  works  against 
it.  Second,  capping  forms  armchair-like  sites  that  initiate  zipper-type  growth  of  new  layers. 
Finally,  capping  permanently  embeds  five-member  rings  into  the  growing  graphene  structure 
and  hence  causes  it  to  curve. 

In  the  Whitesides  and  Frenklach  study  [10],  this  reaction  was  assumed  to  have  the  rate 
coefficient  of  its  zigzag-edge  analog  because  no  rate  data  were  available.  The  importance 
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of  this  reaction  called  for  a  thorough  examination  of  its  kinetics.  To  do  so,  we  employed 
a  combination  of  quantum  chemical  calculations  and  reaction  rate  analysis  to  explore  the 
elementary  steps  of  these  types  of  reactions  and  to  derive  their  rate  coefficients.  To  evaiuate 
the  influence  of  the  embedded  five-member  ring  on  six-member  ring  cyclization,  we  compared 
the  energetics  and  reaction  rates  of  capping  an  embedded  five-member  ring  with  those  of 
capping  a  “regular”  armchair  graphene-edge  site. 

The  rate  coefficient  caiculated  by  us  for  capping  the  five-member-ring  bay  turned  out  to 
be  of  the  same  order  of  magnitude  as  its  prior  estimate,  yet  with  a  substantialiy  different 
activation  energy.  However,  the  calculated  rate  coefficient  for  capping  the  six-member-ring 
bay  is  one  order  of  magnitude  lower  than  that  used  in  prior  KMC  studies. 

To  examine  the  influence  of  the  new  kinetic  data,  we  performed  detailed  kinetic  Monte 
Carlo  simulations  changing  the  rate  coefficients  of  the  bay-capping  reactions  of  the  initiai 
modei.  With  the  new  rate  coefficients,  the  overali  growth  rates  of  graphene  edge  were  iower. 
The  most  dramatic  effect  was  obtained  at  2000  K.  At  the  same  time,  the  change  in  the  rate 
coefficients  did  not  affect  the  phenomenon  of  curving  obtained  for  graphene-edge  evolution 
in  the  previous  study  [10],  as  shown  in  Fig.  2. 


(C) 


Figure  2:  Representative  structures  formed  at  a  simulation  time  of  5  ms  from  coronene 
growth  at  (a)  1500  K,  (b)  2000  K,  and  (c)  2500  K 


This  work  was  presented  at  the  33rd  Symposium  on  Combustion  [74], 

4.1.2  Thermodynamic  stability  of  oxyradicals 

The  oxidative  destruction  and  resistance  to  oxidation  of  soot  are  of  practical  interest,  as  these 
are  primary  pathways  through  which  soot  production  is  alleviated.  However,  experimental 
data  are  scarce.  When  such  data  are  available,  they  typically  have  been  collected  on  a 
composition  of  soot  which  is  not  well  determined  and  dependent  on  many  experimental 
factors,  making  extrapolation  to  new  systems  challenging.  Therefore,  theoretical  methods 
are  ideal  as  they  provide  control  over  modeling  conditions  and  facilitate  extension  to  larger 
and  more  complex  systems. 
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Our  first  objective  is  investigation  of  the  thermodynamic  stability  of  graphene  oxyradi- 
cals,  which  showed  a  systematic  trend  with  position  of  the  chemisorbed  oxygen  atom  on  the 
graphene  edge.  Our  second  objective  is  to  explain  the  observed  trend,  and,  to  this  end,  we 
employ  the  theory  of  aromaticity.  We  chose  the  linear  pentacene  polyaromatic  hydrocarbon 
(PAH)  as  the  model  of  the  graphene  edge.  Because  the  present  focus  is  on  processes  that 
involve  only  the  exterior  of  a  graphene  sheet,  this  model  should  be  adequate.  It  contains  the 
essential  structural  element  involved  in  the  key  processes,  i.e.  the  zigzag  edge  which  can  be 
systematically  extended  by  both  increasing  the  length  of  the  polycyclic  chain  and  by  adding 
extra  layers.  We  undertake  a  systematic  investigation  into  the  thermodynamic  stability  of 
pentacene  oxyradicals  as  a  model  for  the  critical  intermediate  in  graphene  oxidation.  The  re¬ 
sult  of  this  investigation  is  an  intuitive  chemical  model  based  on  aromaticity  and  localization 
that  explains  the  observed  trends. 

We  started  with  a  series  of  linear  polyaromatics,  investigating  the  relative  thermodynamic 
stability  of  pentacene  oxyradicals,  shown  in  Fig.  3.  Oxyradicals  with  O  bonded  to  interior 
rings  were  found  to  be  more  stable  than  those  with  O  bonded  to  exterior  rings  in  the  range 
of  combustion  temperatures.  To  assess  thermodynamic  stability,  we  calculated  the  standard 
Gibbs  free  energies  of  the  four  pentacene  oxyradicals  relative  to  oxyradical  II  with  results 
shown  in  Fig.  4. 

We  then  expanded  the  investigation  to  larger,  two-dimensional  aromatics,  those  shown 
in  Fig.  5.  As  in  the  case  of  the  linear  aromatics,  the  energy  and  hence  the  thermodynamic 
stability  of  larger,  two-dimensional  aromatic  oxyradicals  are  found  to  depend  critically  on 
the  position  of  oxygen  in  the  PAH  edge.  This  is  demonstrated  in  Fig.  6  where  we  plot  the 
standard  enthalpies  of  formation  of  oxyradical  versus  values  of  HOMA  averaged  per  ring  and 
the  standard  deviation  of  individual  HOMA  values.  As  can  be  seen  from  this  figure,  there  is 
a  nearly  linear  correlation  of  the  data,  which  indicates  that  the  radical  stability  depends  on 
the  value  of  HOMA  as  well  as  a  distribution  of  HOMA  over  the  molecular  substrate.  These 
results  indicate  that  oxyradical  stability  correlates  with  local  aromatic  character  of  six-atom 
rings  characterized  by  the  harmonic  oscillator  measure  of  aromaticity  (HOMA)  and  with  the 
distribution  of  HOMA  values  in  molecules.  It  is  demonstrated  that  oxidation  at  the  edge 
has  a  non-local  effect  on  the  structure  of  PAHs  and  leads  to  distinguishable  types  of  HOMA 
patterns  that  are  common  for  both  families  of  PAHs. 

The  interplay  between  local  aromaticity  of  benzene-like  rings  of  graphene  substrates  and 
the  global  aromaticity  of  the  latter  has  led  to  a  simple  explanation  of  the  trends  in  relative 
stability  of  two  families  of  linear  and  two-dimensional  PAH  oxyradicals.  For  any  substrate, 
chemisorption  of  an  oxygen  atom  leads  to  the  loss  of  aromaticity  by  the  corresponding  ring 
and  eliminates  its  contribution  to  inter-ring  conjugation.  As  a  result,  the  relative  stability 
of  linear  oxyradicals  is  controlled  by  fragmentation  of  the  globally  delocalized  7r-electron 
system.  It  mainly  depends  on  the  nature  and  the  amount  of  locally  7r-aromatic  fragments 
formed.  Relative  energies  of  linear  oxyradicals  show  linear  dependency  of  the  cumulative 
HOMA  aromaticity  measure.  This  linear  trend  is  essentially  preserved  in  oxyradicals  of  two- 
dimensional  substrates.  No  fragmentation  but  rather  rearrangement  of  the  pattern  of  local 
aromaticity  occurs  upon  oxidation  in  the  edge.  It  appears  to  lead  to  three  classes  of  local 
aromaticity  patterns.  The  first  class  is  checkerboard-like.  It  is  related  to  the  Clar  structure 
of  coronene  with  highly  aromatic  disjoint  rings.  The  second  class  resembles  a  resonance  of 
Clar  structures  leading  to  cycles  of  six  weakly  aromatic  rings,  hence  the  designations  “Clar 
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Figure  3:  Left  panel:  structures  of  the  pentacene  molecule  (I)  and  pentacene  oxyradicals 
with  the  oxygen  atom  in  different  positions  (II  V).  Connections  between  atoms  are  drawn 
on  the  basis  of  the  interatomic  distances.  Upper  case  Roman  numerals  designate  oxyradicals 
with  different  O  atom  positions;  lower  case  italic  Roman  numerals  designate  six-atom  rings; 
Arabic  numerals  enumerate  C  atoms.  Right  panel:  schematic  representation  of  chemical 
bonding  in  systems  I-V. 
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Cumulative  HOMA 

Figure  4:  ZPE-corrected  relative  energies  of  oxyradicals  II-V  (at  the  B3LYP/6-311G(d,p) 
level  of  theory)  plotted  against  cumulative  HOMA.  The  straight  line  reflects  nearly  linear 
dependence. 


coronene”  and  “superaromatic  coronene” .  Also,  there  is  an  “intermittent”  class  which  can  be 
seen  as  a  mixture  of  the  previous  two.  It  remains  to  be  seen  if  these  three  types  of  patterns 
are  of  general  relevance  and  emerge  in  other  graphene-based  systems.  Their  formation  does 
not  seem  to  be  affected  by  the  nature  of  the  substrate  edges  in  the  present  research,  so  the 
assumption  of  generality  is  plausible.  This  shows  that  a  very  local  event,  such  as  oxidation 
of  the  edge,  has  a  strongly  non-local  effect  on  the  entire  framework  of  C-C  bonds. 

These  results  are  reported  in  Refs.  75,  76,  and  77. 

4.1.3  Kinetics  of  oxyradical  decomposition 

Progress  has  been  made  in  understanding  reaction  mechanisms  of  PAH  formation  and  PAH- 
edge  growth  at  high  temperatures  [4-6,  10].  Oxidation  of  PAH  edges,  while  part  of  the 
overall  growth  processes  in  oxygen-containing  environments,  has  received  substantially  less 
attention.  The  goal  of  the  present  study  is  to  explore  PAH-edge  degradation  mechanisms, 
focusing  on  the  fate  of  PAH  oxyradicals. 

A  number  of  studies  investigated  the  oxidation  of  small  aromatic  molecules  and  oxygen 
chemisorptions  at  selective  sites  of  two-  and  three- ring  aromatics  (for  references  see  Ref.  78). 
According  to  these  studies,  at  high  temperatures  the  reaction  between  a  phenyl  radical  and 
molecular  oxygen  leads  to  the  formation  of  phenoxy  radical,  which  subsequently  decomposes 
to  generate  cyclopentadienyl  radical,  C5H5,  and  carbon  monoxide,  CO,  shown  in  Fig.  7. 
Similarly,  a  naphthoxy  radical  is  produced  from  a  naphthyl  radical  reacting  with  O2  and 
then  decomposes  to  an  indenyl  radical  and  CO. 

Fewer  studies  have  been  carried  out  on  edge  oxidation  of  larger  PAHs  and  those  pertinent 
to  the  present  work  are  all  theoretical.  Radovic  [79]  examined  the  energetics  and  feasible 
pathways  for  the  oxidation  to  CO2  of  several  five-  to  seven-ring  mostly  peri-condensed  PAHs. 
Celnik  et  al.  [13]  computed  the  rate  of  O2  reaction  with  a  PAH  radical  site  generating  an 
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Figure  5:  Structures  of  two-dimensional  PAH  substrates  used  to  form  oxyradicals.  Roman 
numerals  label  position  of  oxygen  in  corresponding  oxyradicals.  Arabic  numerals  label  six- 
atomic  rings. 
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Figure  6:  Standard  enthalpy  of  formation,  in  kcal/mol,  of  two-dimensional  oxyradicals  versus 
average  values  of  HOMA  and  its  standard  deviation,  ohoma-  Also  shown  is  a  plane  that  is 
a  fit  through  the  data  points.  The  inset  displays  the  same  figure  but  rotated  to  show  the 
projection  normal  to  the  plane. 
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Figure  7:  Phenyl  oxidation  by  O2. 
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oxyradical  and  an  oxygen  atom  for  benzene,  naphthalene,  anthracene,  phenanthrene,  pyrene, 
and  benzo [a] anthracene  at  the  B3LYP/6-31G(d)  level  of  theory.  Their  results,  in  accord 
with  an  earlier  study  of  Kunioshi  et  al.  [80] ,  showed  that  the  oxidation  rates  are  similar  for 
all  radical  sites  studied  except  for  armchair  sites,  which  exhibited  higher  activation  energy. 
However,  the  decomposition  of  the  oxyradicals  generated  was  not  investigated  in  their  study. 

All  of  these  studies  suggest  that  oxyradicals  are  not  only  the  major  intermediates  during 
oxidation  of  small  aromatics,  but  also  of  large  PAHs  and  graphene  edges.  It  is  of  inter¬ 
est,  therefore,  to  investigate  the  fate  of  PAH  oxyradicals  and  shed  light  on  the  evolution 
of  graphene  edges  and  soot  surfaces  during  oxidation.  Specifically,  we  want  to  determine 
whether  graphene-edge  oxyradicals  decompose  similarly  to  phenoxy  and  naphthoxy  radicals 
and,  if  so,  with  what  rates.  We  answer  these  questions  by  performing  a  theoretical  study, 
following  the  previous  work  on  the  thermodynamic  stability  of  graphene  oxyradicals,  at  the 
same  level  of  theory. 

First,  we  focus  on  zigzag  edges.  For  this  we  select  pentacene,  as  in  our  prior  thermo¬ 
dynamic  study.  We  examine  the  reaction  kinetics  of  four  pentacene  oxyradicals,  shown  in 
Fig.  8,  as  well  as  that  of  phenoxy  radical,  by  solving  energy-transfer  master  equations. 

We  examined  four  pentacene  oxyradicals  distinguished  by  O-atom  bonding  to  carbon 
atoms,  denoted  I,  II,  III,  and  IV,  respectively;  these  structures  are  displayed  in  Fig.  8. 
Although  the  four  pentacene  oxyradicals  look  similar  in  structure,  they  are  quite  different  in 
potential  energy  and  thermodynamic  stability.  Pentacene  oxyradicals  I,  II,  and  III  are  17.1, 
14.1,  and  2.9  kcal/mol  higher  in  potential  energy,  respectively,  than  oxyradical  IV  at  the 
B3LYP/6-311G(d,p)  level.  We  showed  earlier  [75]  that  below  1000  K  the  thermodynamic 
stability  of  the  four  oxyradicals  follows  the  trend  of  their  relative  potential  energies,  i.e. 
the  order  IV  >  III  >  II  >  I,  while  above  1000  K  oxyradical  III  becomes  thermodynamically 
more  stable  than  IV  because  of  the  larger  entropy  contribution  to  the  Gibbs  free  energy.  The 
relative  energies  and  stabilities  can  be  explained,  as  for  the  thermodynamics  case  Fig.  8,  by 
the  different  fragmentation  of  the  delocalized  7r-electron  system  of  the  precursor  pentacene 
molecule  arising  from  the  different  locations  of  the  chemisorbed  oxygen  atom. 

Examination  of  the  decomposition  path  of  the  phenoxy  radical  indicates  that  the  oxygen 
atom,  01,  and  the  three  carbon  atoms,  C2,  C3,  and  C4  (labeled  in  Fig.  8)  are  involved  in 
breaking  existing  and  forming  new  C-C  bonds.  It  is  worthwhile  comparing  bond  lengths 
of  these  four  atoms  for  pentacene  oxyradicals  with  those  for  phenoxy  radical.  As  shown  in 
Fig.  8,  the  bond  lengths  or  distances  between  Ol  and  C2,  C2  and  C3,  C2  and  C4,  and  C3 
and  C4  of  oxyradical  I  are  closest  to  those  of  the  phenoxy  radical  among  the  four  pentacene 
oxyradicals.  For  instance,  the  distance  between  C3  and  C4  is  in  the  order  of  phenoxy  <  I 
<  II  <  III  <  IV.  Therefore,  the  similarity  in  structure  of  the  pentacene  oxyradicals  to  the 
phenoxy  radical  is  in  the  order  of  I  >  II  >  III  >  IV.  An  oxyradical  with  a  shorter  distance 
between  C3  and  C4  is  more  likely  to  decompose  in  the  same  way  as  the  phenoxy  radical  and 
to  generate  a  cyclic  intermediate. 

We  explored  the  potential  energy  surfaces  at  the  B3LYP/6-311G(d,p)  level  for  the  ther¬ 
mal  decomposition  of  the  four  pentacene  oxyradicals.  To  facilitate  analysis  of  the  pentacene 
oxyradical  systems,  we  also  computed  the  phenoxy  reaction  system  at  the  same  level  of 
theory.  A  detailed  report  of  the  computed  molecular  properties — electronic  energies,  geome¬ 
tries,  and  vibrational  frequencies — for  the  species  involved  is  given  in  Ref.  78.  Our  results 
indicate  that  the  decomposition  pathways  of  pentacene  oxyradicals  I  and  II  leading  to  CO 
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Figure  8:  Structures  of  phenoxy  and  pentacene  oxyradicals 
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generation  are  very  similar  to  that  of  phenoxy  radical,  while  CO  is  unlikely  to  be  produced 
following  the  same  pathways  for  oxyradicals  III  and  IV. 


TS5-3 


Figure  9:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  phenoxy  radical 
at  the  B3LYP/6-311G(d,p)  level.  Energies  (in  kcal/mol)  at  0  K  relative  to  reactants  are 
shown  for  each  transition  state  and  local  minimum. 


Figure  9  shows  the  minimum  potential-energy  paths  for  the  decomposition  of  the  phenoxy 
radical  following  different  pathways.  One  proceeds  through  an  electro-cyclic  mechanism, 
forming  the  cyclic  intermediate  2,  followed  by  C-C  bond  cleavage  to  produce  intermediate  3. 
Elimination  of  CO  from  3  leads  to  the  main  products,  4  and  CO.  The  other  pathway  proceeds 
through  the  ring  opening,  i.e.,  by  breaking  either  the  C2-C3  or  C2-C4  bond.  As  noted,  the 
ring  opening  requires  much  higher  activation  energy  than  the  cyclization  mechanism,  in 
agreement  with  earlier  studies  (for  references  see  [78]). 

As  shown  in  Figs.  10  and  11,  pentacene  oxyradicals  I  and  II  follow  pathways  similar 
to  phenoxy  radical.  However,  the  reactions  generating  products  4  and  CO  from  pentacene 
oxyradicals  I  and  II  are  less  endothermic,  by  13.9  and  11.0  kcal/mol,  respectively,  than  those 
of  phenoxy  radical.  Compared  with  phenoxy  radical,  pentacene  oxyradical  I  has  a  lower 
barrier  (TS1-2)  to  generate  intermediate  2,  while  oxyradical  II  has  a  higher  barrier.  The 
TS1-2  energy  of  the  pentacene  oxyradical  I  system  is  lower  than  that  of  II  by  10.4  kcal/mol, 
which  contributes  to  a  faster  decomposition  rate  of  oxyradical  I. 

Similar  to  phenoxy  radical,  pentacene  oxyradicals  I  and  II  each  have  distinct  possible 
ring-opening  pathways  by  breaking  either  of  the  two  neighboring  C-C  bonds.  For  phenoxy 
radical,  these  two  pathways  are  the  same.  However,  for  pentacene  oxyradicals  I  and  II, 
the  two  possible  ring-opening  pathways  are  different.  One  of  these  pathways,  via  TS1-5,  is 
comparable  in  energy  to  that  of  phenoxy  radical.  The  other  one,  via  TS1-6,  is  much  higher 
in  energy,  which  makes  the  fate  of  6  essentially  unimportant. 

In  contrast  to  pentacene  oxyradicals  I  and  II,  no  cyclic  intermediates  were  found  for 
pentacene  oxyradicals  III  and  IV.  A  potential  energy  scan,  in  which  the  C3  to  C4  distance 
was  varied  and  other  geometrical  parameters  allowed  to  relax,  was  performed  to  explore 
possible  reaction  pathways  for  III  and  IV  along  the  reaction  coordinate  to  form  a  cyclic 
intermediate,  as  shown  in  Fig.  12.  For  oxyradical  III,  as  the  two  interacting  carbon  atoms 
come  close  to  form  a  bond,  only  a  local  inflection  point  was  found,  indicating  that  a  local 
minimum  either  does  not  exist  or  the  potential  energy  well  is  too  shallow.  In  the  case  of 
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Figure  10:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  pentacene 
oxyradical  I  at  the  B3LYP/6-311G(d,p)  level.  Energies  (in  kcal/mol)  at  0  K  relative  to 
reactants  are  shown  for  each  transition  state  and  local  minimum. 


Figure  11:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  pentacene 
oxyradical  II  at  the  B3LYP/6-311G(d,p)  level.  Energies  (in  kcal/mol)  at  0  K  relative  to 
reactants  are  shown  for  each  transition  state  and  local  minimum. 
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Figure  12:  Relaxed  potential  energy  scan  for  pentacene  oxyradicals  III  (left)  and  IV  (right) 
at  the  B3LYP/6-311G(d,p)  level  of  theory. 


Figure  13:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  pentacene 
oxyradical  III  (left)  and  IV  (right)  obtained  by  breaking  the  C2-C3  and  C2-C4  bonds  adjacent 
to  the  C=0  bonds,  respectively,  at  the  B3LYP/6-311G(d,p)  level.  Energies  (in  kcal/mol)  at 
0  K  relative  to  reactants  are  shfown  for  each  transition  state  and  local  minimum. 


31 


oxyradical  IV,  the  potential  energy  continues  to  increase  as  the  two  carbon  atoms  get  closer. 
Fig.  13  depicts  the  decomposition  of  oxyradicals  III  and  IV  through  their  respective  ring¬ 
opening  pathways.  In  both  of  these  cases,  intermediates  5  and  6  have  very  high  potential 
energies  relative  to  reactants  with  potential-energy  barriers  over  100  kcal/mol.  Therefore,  it 
is  unlikely  for  oxyradicals  III  and  IV  to  produce  CO  through  reaction  pathways  similar  to 
those  of  oxyradicals  I  and  II.  Furthermore,  the  ring-opening  pathways  of  oxyradicals  III  and 
IV  involve  internal  torsional  rotations  and  therefore  form  non-planar  molecular  structures, 
which  is  unlikely  to  happen  for  multi-layer  PAHs. 

Based  on  the  potential  energy  surfaces  and  the  molecular  properties  calculated  at  the 
B3LYP/6-311G(d,p)  level,  we  computed  the  rate  coefficients  of  the  thermal  decomposition  of 
phenoxy  radical  and  the  pentacene  oxyradicals.  For  phenoxy  and  pentacene  oxyradicals  I  and 
II,  we  found  that  the  main  decomposition  products  were  carbon  monoxide  and  species  4,  and 
the  production  of  other  stable  species  was  negligible  for  the  conditions  studied.  Therefore, 
we  may  write  the  reactions  as  shown  in  Fig.  14.  For  pentacene  oxyradicals  I  and  II,  the 


Figure  14:  Overall  decomposition  reactions  of  aromatic  oxyradicals. 


ring-opening  pathway  through  TS1-6  can  be  neglected  as  its  rate  coefficient  is  orders  of 
magnitude  smaller  than  that  through  TS1-5.  Assuming  intermediate  species  2,  3,  and  5  are 
in  steady  state,  the  high-pressure- limits  of  these  rate  coefficients,  k^ ,  can  be  calculated  for 
the  overall  reactions  Rl,  R2,  and  R3, 

khoo=  1.84  x  107TL93  exp  (-25819/T)  s"1 
k2, oo  =  9.23  x  107TL61  exp  (-25078 /T)  s"1 
h)OQ=  1.11  x  107TL97  exp  (-29070/T)  s”1 

Analysis  of  the  computed  reaction  rates  indicated  that  the  cyclic  pathway  and  the  ring¬ 
opening  pathway  through  TS1-5  both  contribute  to  ,  and  that  the  rate  coefficients  of 
steps  1— >2  and  1— >5  are  smaller  than  those  of  the  subsequent  steps. 

The  rate  coefficients  at  finite  pressures  were  computed  by  solving  the  master  equations 
using  the  Multiwell  suit  of  codes,  and  the  results  are  displayed  in  Figs.  15,  16,  and  17. 

To  compare  our  results  with  those  in  the  literature  on  phenoxy  decomposition,  Fig.  15  dis¬ 
plays  two  sets  of  reported  experimental  results  and  our  computed  k\  as  a  function  of  tempera¬ 
ture  and  pressure.  One  experimental  study  is  that  of  Lin  and  Lin  [81] ,  who  investigated  Rl  at 
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Figure  15:  Comparison  of  computed  k\  with  experimental  rate  coefficients. 


Figure  16:  Rate  coefficients  of  CO  formation  from  pentacene  oxyradicals  I  and  II  as  a  function 
of  temperature  and  pressure.  Solid  lines:  pentacene  oxyradical  I;  dashed  fines:  pentacene 
oxy radical  II. 
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Figure  17:  Rate  coefficients  of  CO  formation  from  phenoxy  and  pentacene  oxyradicals  I  and 
II.  Solid  fines:  the  high-pressure  limit;  dashed  fines:  1  atm. 


1000-1580  K  and  0.4-0. 9  atm  using  anisole  (C6H5OCH3)  and  allyl  phenyl  ether  (C6H5OCH2) 
as  precursors.  These  authors  reported  a  rate  coefficient  of  10n'4=to  :2exp(— 22100  ±  450/T) 
s_1.  The  other  experimental  result  comes  from  a  shock-tube  study  of  phenyl  radicals  reacting 
with  molecular  and  atomic  oxygen  by  Frank  et  al.  [82],  who  determined  the  rate  coefficient 
to  be  7.4  x  10nexp(— 22070/T)  s_1  in  the  range  of  1020-1190  K  and  1.3-2. 5  bar. 

Compared  with  the  experimental  studies,  our  computed  results  are  within  the  uncertainty 
ranges  of  the  measured  values.  This  agreement  can  serve  as  validation  of  the  present  method 
and  procedures.  We  note  that  the  computed  k\  is  both  temperature  and  pressure  dependent. 
As  expected,  k\  increases  as  temperature  or  pressure  rises.  At  relatively  low  temperatures, 
k\  is  close  to  the  high-pressure  limit,  A-ioo,  when  the  pressure  is  larger  than  1  atm.  At 
intermediate  temperatures,  it  is  in  the  falloff  region.  At  high  temperatures,  it  approaches 
the  low-pressure  limit  due  to  high  internal  energy  at  these  temperatures.  The  pressure 
dependence  of  k\  may  explain  the  discrepancy  between  the  experimental  studies,  since  the 
measurements  of  Frank  et  al.  [82]  were  carried  out  at  higher  pressures  than  those  of  Lin  and 
Lin  [81]. 

Figure  16  shows  the  computed  rate  coefficients  A;2  and  /»;3  of  CO  formation  from  the 
decomposition  of  pentacene  oxyradicals  I  and  II,  respectively.  Similar  to  phenoxy  radical, 
the  deviation  of  the  rate  coefficients  from  their  respective  high-pressure  limits  increases  with 
temperature.  We  also  observe  that  in  the  temperature  range  of  1500-2500  K  and  the  pressure 
range  of  0.01-10  atm,  the  computed  values  of  k2  are  1.6  to  8.2  times  those  of 

Figure  17  displays  the  rate  coefficients  of  CO  formation  from  the  decomposition  of  phe¬ 
noxy  radical  and  pentacene  oxyradicals  I  and  II  computed  at  the  high-pressure  limit  and  at 
1  atm.  As  can  be  seen  in  Fig.  17,  the  rate  coefficients  are  pressure  dependent  for  all  three 
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species.  At  the  high-pressure  limit,  k\ ^  is  larger  than  both  k2oo  and  hioo ,  while  at  1  atm, 
h  is  smaller  than  k2  and  close  to  h> .  The  pressure  dependence  is  weaker  for  the  pentacene 
oxyradicals  than  for  the  phenoxy  radical,  as  expected  for  large  molecules. 

The  results  obtained  indicate  that  the  decomposition  rate  is  determined  by  the  location 
of  the  oxygen  atom  in  the  pentacene  oxyradicals.  Oxyradicals  with  oxygen  attached  to  inner 
rings  are  kinetically  more  stable  than  those  with  oxygen  attached  to  outer  rings.  The  latter 
can  generate  carbon  monoxide  at  rates  comparable  to  phenoxy  radicals,  while  CO  is  unlikely 
to  be  produced  from  oxyradicals  with  oxygen  bonded  to  inner  rings. 

The  kinetics  results  for  the  zigzag  edges  are  published  in  [78]. 

Finally,  we  computed  decomposition  rates  of  graphene  oxyradicals  whose  oxygen  atom  is 
located  at  armchair  sites,  compare  these  rates  to  those  at  zigzag  sites  dicussed  above.  We 
then  show  that  kinetic  stability  of  oxyradicals,  similar  to  their  thermodynamic  stability,  is 
correlated  with  substrate  aromaticity. 

To  study  the  decomposition  of  armchair  edge  oxyradicals,  five  molecules  were  selected; 
these  structures  are  shown  in  Fig.  18.  Phenanthrene  oxyradicals  I  and  II  were  selected 


Specie  name 


Phenanthrene  I 


Phenanthrene  II 


Phenanthrene  III 


Benzoperylene 

Extended- 

armchair 


Structure 


Figure  18:  Armchair  oxyradical  structures. 


to  examine  the  differences  in  decomposition  rate  between  the  two  principal  armchair  sites. 
Phenanthrene  oxyradical  III  was  included  to  provide  comparison  with  the  previously  studied 
zigzag  site  decomposition  [78] .  Benzoperylene  oxyradical  and  the  extended-armchair  oxyrad¬ 
ical  have  oxygen  located  on  an  armchair  site  that  is  the  same  as  in  phenanthrene  II,  but  part 
of  a  larger  substrate.  The  intent  of  including  these  larger  substrates  is  to  examine  if  and  by 
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how  much  the  decomposition  rate  for  the  same  edge  site  changes  with  substrate  size. 

Phenanthrene  oxyradicals  I,  II,  and  III  decompose  via  similar  pathways  and  intermediates 
(Figs.  19,  20,  and  21). 


Figure  19:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  phenanthrene 
oxyradical  I  at  the  B3LYP/6-311G(d,p)  level.  Energies  are  in  kcal/mol  at  0  K  relative  to 
the  reactant. 


The  three  pathways  contain  a  cyclic  intermediate  2,  and  all  have  at  least  one  ring  open¬ 
ing  intermediate  5.  For  phenanthrene  oxyradical  I,  two  five-member  rings  with  attached 
CO,  intermediates  3  and  7,  were  found.  For  phenanthrene  oxyradicals  II  and  III,  only  one 
such  intermediate,  3,  was  discovered.  An  additional  channel  was  found  for  phenanthrene 
oxyradical  II,  namely  TS  2-4,  which  connects  intermediate  2  directly  to  product  4. 

These  pathways  are  similar  to  those  computed  for  the  decomposition  of  pentacene  zigzag 
corner  oxyradicals,  which  included  the  cyclic  (2),  attached-CO  (3),  and  ring-opening  (5) 
intermediates.  Furthermore,  the  products  of  decomposition  are  also  comparable  in  that 
all  contain  a  five-member  ring.  Barrier  energies  of  the  corner-pentacene  and  phenanthrene 
oxyradical  decomposition  pathways  are  also  very  close  to  each  other,  often  to  within  3- 
4  kcal/mol.  Such  parallels  in  the  reaction  pathways  indicate  that  the  armchair-edge  and 
corner  zigzag-edge  decompositions  are  comparable  and  could  be  classified  as  equivalent  types 
of  oxidation  reactions. 

Phenanthrene  oxyradicals  I,  II,  and  III  are  all  single  row  structures.  To  test  the  influ¬ 
ence  of  substrate  size  on  decomposition  kinetics,  the  benzoperylene  and  extended-armchair 
oxyradicals  were  studied;  both  have  oxygen  on  a  site  similar  to  that  in  phenanthrene  II. 
The  phenanthrene  II  and  benzoperylene  decomposition  pathways,  shown  in  Figs.  20  and 
22,  respectively,  contain  the  same  intermediates,  but  the  energy  of  the  highest  transition 
state  barrier  along  the  lowest  potential  energy  path  (TS  1-2)  is  larger  by  7  kcal/mol  for 
benzoperylene.  As  shown  below,  this  leads  to  an  almost  tenfold  slower  decomposition  rate 
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Figure  20:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  phenanthrene 
oxyradical  II  at  the  B3LYP/6-311G(d,p)  level.  Energies  are  in  kcal/mol  at  0  K  relative  to 
the  reactant. 


Figure  21:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  phenanthrene 
oxyradical  III  at  the  B3LYP/6-311G(d,p)  level.  Energies  are  in  kcal/mol  at  0  K  relative  to 
the  reactant. 
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for  benzoperylene  oxyradical  relative  to  phenanthrene  II  oxyradical.  The  extended-armchair 
decomposition  pathway  differs  from  the  phenanthrene  II  and  benzoperylene  pathways.  The 
extended-armchair  pathway,  shown  in  Fig.  23,  does  not  contain  the  cyclic  intermediate  2, 
but  it  does  have  the  ring-opening  intermediate  5.  In  this  case,  rather  than  proceeding  from 
1  to  2  to  3  as  with  the  other  structures,  the  extended- armchair  oxyradical  goes  directly  from 
1  to  3.  The  difference  in  potential  energy  between  the  highest  barrier  on  the  lowest  energy 
path  of  the  extended-armchair  (TS  1-3)  and  that  of  benzoperylene  (TS  1-2)  is  2.8  kcal/mol. 
This  difference  is  counteracted  by  a  11.7  kcal/mol  higher  ring-opening  barrier  (TS  1-5)  of 
benzoperylene  compared  to  that  of  the  extended- armchair.  These  opposing  trends  nearly 
cancel  and,  as  shown  in  the  next  section,  the  extended-armchair  and  benzoperylene  oxyrad- 
icals  decompose  at  similar  rates.  Using  the  above  results,  reaction  rates  were  calculated  for 


Figure  22:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  benzoperylene 
oxyradical  at  the  B3LYP/6-31lG(d,p)  level.  Energies  are  in  kcal/mol  at  0  K  relative  to  the 
reactant. 


the  five  reactions,  shown  in  Fig.  24. 

The  principal  products  of  oxyradical  decompositions  were  the  embedded  five-member  ring 
and  CO.  The  accumulation  of  other  species  was  negligible  for  the  conditions  studied.  The 
rate  coefficients  at  specific  values  of  pressure  and  temperature  were  computed  by  solving  the 
master  equations  using  the  Multiwell  suite  of  codes.  All  of  these  reactions  displayed  similar 
pressure  and  temperature  dependence. 

In  observing  the  pressure  dependence,  it  was  desirable  to  calculate  the  high-pressure- 
limit  rate  coefficients  for  each  of  the  overall  reactions  Rla-R5a.  To  accomplish  this  aim,  we 
assumed  that  all  intermediates  of  an  individual  reaction  system  are  in  steady  state.  The  fitted 
Arrhenius  expressions  for  the  high-pressure-limit  rate  coefficients  are  reported  in  Fig.  25. 

Inspection  of  rate  coefficients  computed  for  the  three  phenanthrene  oxyradicals,  displayed 
in  Fig.  26,  indicates  that  phenanthrene  I,  with  oxygen  furthest  to  the  outside,  decomposes 


Figure  23:  Minimum  potential  energy  paths  for  the  thermal  decomposition  of  the  extended- 
armchair  oxyradical  at  the  B3LYP/6-311G(d,p)  level.  Energies  are  in  kcal/mol  at  0  K 
relative  to  the  reactant. 


faster  than  phenanthrene  II,  whereas  the  latter  decomposes  about  as  fast  as  phenanthrene 
III.  The  correlation  of  these  rates  with  the  molecular  structure  of  the  three  phenanthrene 
oxyradicals  can  be  rationalized  by  counting  “free  edges”,  a  concept  introduced  in  Ref.  78. 
A  “free  edge”  classifies  the  connectivity  of  the  C-0  site  in  the  context  of  the  multi-ring 
aromatic  structure.  For  example,  in  phenanthrene  I,  C-0  is  connected,  on  both  sides,  to 
carbon  atoms  that  are  members  of  only  one  ring;  we  classify  each  of  these  edges  as  free 
edges.  On  the  other  hand,  C-0  in  phenanthrenes  II  and  III  is  connected  on  one  side  to  a  C 
atom  that  is  a  member  of  two  rings,  and  on  the  other  side  to  a  C  atom  that  is  a  member 
of  one  ring.  Therefore,  phenanthrenes  II  and  III  contain  only  one  free  edge,  with  the  non- 
free  edge  connected  to  the  carbon  atom  that  is  a  member  of  two  rings.  Using  this  simple 
structural  classification  enables  us  to  rationalize  the  relative  ordering  of  the  decomposition 
rates  among  phenanthrene  oxyradicals  I,  II,  and  III.  Phenanthrene  I  has  two  free  edges  and 
decomposes  fastest,  whereas  II  and  III,  each  having  just  one  free  edge,  decompose  slower 
than  I  and  at  nearly  comparable  rates. 

In  Fig.  27  we  compare  the  decomposition  rates  of  phenanthrene  II,  benzoperylene,  and 
the  extended-armchair  oxyradicals.  Inspection  of  these  Endings  indicates  that  the  increase 
in  substrate  size,  from  phenanthrene  to  benzoperylene,  substantially  reduces  the  computed 
reaction  rate,  whereas  the  rate  seems  to  level  off  with  further  increase  in  size,  from  benzop¬ 
erylene  to  extended- armchair.  The  benzoperylene  oxyradical  decomposes  nearly  an  order  of 
magnitude  slower  than  the  phenanthrene  II  oxyradical.  This  is  primarily  due  to  the  higher 
potential  energy  of  TS1-2  in  the  benzoperylene  pathway  compared  to  that  of  phenanthrene 
II.  Comparing  the  extended- armchair  energetics  to  that  of  benzoperylene  indicates  that  al¬ 
though  the  lowest  energy  path  for  the  extended-armchair  has  a  higher  barrier  (TS  1-3)  than 
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Figure  24:  Overall  decomposition  reactions  of  aromatic  armchair  oxyradicals. 
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Reaction 

Rate  coefficient  (s'1) 

Rla 

1.08xl02r3-28exp(— 21866/7) 

R2a 

4.00xl0_274,38exp(— 24624/7) 

R3a 

6.42xl0lo7looexp(— 29183/7) 

R4a 

3.78xl0472-59exp(— 30151/7) 

R5a 

4.47xl0127O39exp(— 33796/7) 

Figure  25:  High-pressure-limit  rate  coefficients. 


benzoperylene  (TS  1-2),  it  has  a  substantially  lower  barrier  to  ring  opening,  thus  making 
its  decomposition  rate  competitive  with  that  of  benzoperylene  oxyradical  decomposition. 
These  results  may  indicate  that  a  two-row  graphene  structures  may  be  sufficiently  large  to 
represent  an  edge  reaction  for  an  arbitrarily  large  graphene  substrate. 

Comparing  the  decomposition  rates  of  the  armchair-edge  and  zigzag-edge  oxyradicals 
highlights  the  resemblance  between  the  decomposition  rates  of  the  armchair  and  corner 
zigzag  edges.  Figure  28  displays  the  decomposition  rates  of  phenanthrene  I  and  II  as  well  as 
those  of  pentacene  I  and  II  from  our  earlier  zigzag  study.  Focusing  first  on  the  1-atm  results 
displayed  in  Fig.  28,  we  see  strong  correlation  between  the  phenanthrene  and  pentacene 
decomposition  rates,  namely  phenanthrene  I  and  pentacene  I  decomposition  rates  are  very 
close  to  one  another,  as  are  the  decomposition  rates  of  phenanthrene  II  and  pentacene  II. 
Structurally,  phenanthrene  I  is  similar  to  pentacene  I — both  have  free  edges  on  each  side  of 
the  C-0  bond.  Likewise,  phenanthrene  II  structurally  resembles  pentacene  II,  with  a  free 
edge  on  only  one  side  of  the  C-0  bond.  Again,  the  structures  with  the  same  number  of  free 
edges  decompose  at  similar  rates,  and  the  structures  with  two  free  edges  decompose  faster 
than  those  with  one  free  edge.  However,  at  the  high-pressure  limit  the  correlation  between 
the  pentacene  and  phenanthrene  oxyradical  decomposition  rates  is  weaker  (Fig.  28).  The 
general  trend  within  the  different  types  of  molecules  still  holds,  as  oxyradical  structures 
with  more  free  edges  decompose  at  a  faster  rate  than  those  with  just  one  free  edge.  This 
“free-edge”  analysis  rationalizes  one  of  the  key  results  of  the  present  study,  namely,  that 
phenanthrene  oxyradicals  and  corner  pentacene  oxyradicals  decompose  at  similar  rates. 
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Figure  26:  Rate  coefficients  of  the  decomposition  of  phenanthrene  oxyradicals  I,  II  and  III. 
The  solid  fines  are  the  high-pressure- limit  values  and  the  dashed  lines  are  the  rate  coefficients 
at  1  atm. 
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Figure  27:  Rate  coefficients  of  the  decomposition  of  phenanthrene  oxyradical  II  (black  fines), 
benzoperylene  oxyradical  (black  lines),  and  the  extended-armchair  oxyradical  (gray  fines). 
The  solid  fines  are  the  high-pressure- limit  values  and  the  dashed  lines  are  the  rate  coefficients 
at  1  atm. 
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Figure  28:  Rate  coefficients  of  the  decomposition  of  phenanthrene  oxyradicals  I  and  II  (black 
lines)  and  pentacene  oxyradicals  I  and  II  (gray  lines).  The  solid  lines  are  the  high-pressure- 
limit  values  and  the  dashed  lines  are  the  rate  coefficients  at  1  atm. 
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We  now  turn  to  examing  the  correlation  of  the  computed  reaction  rates  with  the  molec¬ 
ular  structure  of  the  oxyradicals.  Polyaromatic  hydrocarbons  are  structurally  related  to  the 
benzene  molecule.  This  connection  implies  strong  similarity  in  electronic  structure  that  char¬ 
acterizes  aromatic  systems,  namely,  that  of  delocalized  bonding.  However,  in  the  analysis  of 
chemical  properties  of  large  PAHs  studied  as  models  of  extended  systems,  such  as  soot  and 
graphene,  this  connection  is  rarely  exploited.  It  is  reasonable  to  examine  thermodynamic  and 
kinetic  aspects  of  reactivity  of  large  PAHs  on  the  basis  of  their  local  and  global  aromaticity. 
Such  a  connection  can  lead  to  formulation  of  simple  quantitative  or  semi-quantitative  pre¬ 
dictive  models  that  appeal  to  chemical  intuition  and  help  to  rationalize  structure-property 
relationships  of  the  target  systems — soot  and  graphene. 

Aromaticity  can  be  defined  as  the  extra  stabilization  of  a  system  featuring  delocalized 
bonds  with  respect  to  a  fictitious  system  with  completely  localized  bonds.  In  practice,  as¬ 
sessment  of  the  resonance  energy  is  a  difficult  task,  especially  for  molecules  that  require 
multiple  Kekule  structures  in  the  resonance  description.  Multiple  measures  of  aromaticity 
have  been  proposed  in  order  to  simplify  the  task  and  connect  the  analysis  to  observable 
properties.  Considering  the  computational  challenges  associated  with  theoretical  treatment 
of  large  PAHs  it  appears  reasonable  to  choose  the  simplest  model  of  aromaticity  that  uti¬ 
lizes  the  observation  that  aromatic  bonding  is  strongly  associated  with  equalization  of  bond 
lengths  due  to  formation  of  a  delocalized  bonding  framework.  In  the  Harmonic  Oscillator 
Model  of  Aromaticity  (HOMA)  employed  in  our  earlier  thermodynamics  study,  one  evaluates 
the  mean  deviation  of  bond  lengths  in  a  system  under  consideration  with  respect  to  bond 
lengths  in  benzene.  It  is  further  scaled  such  that  HOMA  =  0  for  the  Kekule  form  of  benzene 
and  HOMA  =  1  for  the  aromatic  form. 

Our  previous  studies  demonstrated  that  a  simple  correlation  exists  between  cumulative 
HOMA  of  oxyradicals  and  their  relative  thermodynamic  stability.  The  origin  of  this  rela¬ 
tionship  is  fragmentation  of  the  delocalized  7r-bonding  framework  due  to  its  disruption  upon 
chemisorption  of  an  O  atom  at  a  graphene  radical  site.  It  is  reasonable  to  assume  that 
the  same  reasoning  applies  to  energetic  characteristics  of  oxyradical  decomposition.  In  the 
present  case,  HOMA  can  be  used  for  reactivity  analysis  for  reaction  kinetics.  Considering 
the  diversity  of  factors  controlling  reaction  kinetics,  for  example,  the  multiple  intermediates 
and  transition  states,  we  narrowed  the  range  of  analyzed  characteristics  down  to  a  single 
one,  namely  the  relative  energy  of  the  first  transition  state.  Barrier  height  is  a  relative  value, 
so  it  is  logical  to  consider  its  relationship  to  a  relative  descriptor  of  aromatic  bonding,  that 
can  be  defined  as  the  difference  between  HOMA  of  the  benzene-like  ring  undergoing  CO 
expulsion  of  the  reactant  and  HOMA  of  the  transition  state,  i.e.,  HOMArei  =  HOMAts  - 
HOMAreactant-  A  simple  correlation  in  the  data  can  be  seen  by  plotting  (Fig.  29)  the  results 
obtained  for  both  zigzag  and  armchair  cases.  The  smallest  barrier  height  is  associated  with 
the  smallest  absolute  HOMA  change,  thus  suggesting  that  CO  expulsion  should  proceed  via 
the  route  that  involves  the  smallest  rearrangement  of  the  conjugated  bonding  framework. 

4.2  DNS  Study  of  Soot  Formation  and  Evolution 

A  DNS  study  was  performed  to  provide  insight  into  the  processes  controlling  soot  forma¬ 
tion  and  growth  in  turbulent  flames,  while  considering  state-of-the-art  chemistry  and  soot 
models  developed  under  this  project.  A  two-dimensional  turbulent  nonpremixed  flame  sub- 
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Figure  29:  First  barrier  heights  of  oxyradical  decompositions  against  relative  HOMA  of 
respective  transition  states  (TS  1-2). 
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ject  to  decaying  isotropic  turbulence  was  selected  to  make  the  simulation  computationally 
tractable.  This  study  distinguishes  it-  self  from  previous  work  [28,  29]  by  considering  finite 
rate  chemistry  from  fuel  oxidation  to  the  formation  of  PAH  and  a  detailed  soot  model  based 
on  elementary  physical  processes  and  rates  rather  than  using  a  semi-empirical  approach. 
Finite  rate  chemistry  allows  for  the  investigation  of  the  effect  of  turbulent  mixing  on  the 
formation  and  yield  of  key  soot  precursors  such  as  acetylene  and  naphthalene. 

The  complete  formulation  of  the  flow  governing  equations  and  simulation  numerical  meth¬ 
ods  used  in  the  DNS  are  provided  in  Ref.  24,  but  are  omitted  here  to  focus  more  directly  on 
the  study’s  implications  for  LES  modeling. 

The  detailed  chemical  mechanism  used  to  generate  the  LES  flamelet  tables  (see  Sec.  3.3) 
must  be  reduced  to  be  feasible  for  finite  rate  chemistry  computations.  A  multi-step  approach 
combining  automatic  and  manual  reduction  techniques  was  used  to  reduce  the  size  of  the 
chemical  mechanism  and  make  it  more  affordable  for  DNS.  Throughout  the  reduction  process, 
emphasis  was  placed  on  accurately  predicting  the  combustion  characteristics  of  both  n- 
heptane  (C7H16)  and  toluene  (CyHg),  two  major  components  typically  found  in  surrogate 
fuel  formulations.  The  final  reduced  mechanism  contains  only  47  species  and  290  reactions. 
It  was  validated  against  a  series  of  laminar  test  cases,  including  laminar  burning  velocities,  a 
rich  premixed  flame,  and  a  diffusion  flame.  Soot  was  modeled  using  the  approaches  described 
in  Sec.  3.2.1. 

4.2.1  Initial  Conditions  and  Simulation  Parameters 

The  two-dimensional  computational  domain  consists  of  a  square  of  size  L.  Periodic  bound¬ 
ary  conditions  are  applied  in  both  the  horizontal  (x)  and  vertical  (y)  directions,  effectively 
resulting  in  a  constant  volume  system.  At  the  onset  of  the  simulation,  a  horizontal  strip 
of  fuel  is  surrounded  by  oxidizer.  The  fuel  stream  consists  of  n-heptane  diluted  with  84.4% 
(by  volume)  nitrogen  at  a  temperature  of  300  K.  The  oxidizer  stream  consists  of  air  (21% 
oxygen  and  79%  nitrogen)  also  at  300  K.  The  significant  nitrogen  dilution  in  the  fuel  stream 
is  necessary  to  reduce  peak  soot  volume  fraction  to  a  level  for  which  radiative  heat  transfer 
from  soot  particles  and  its  effect  on  the  gas-  phase  can  be  safely  neglected.  The  compositions 
of  fuel  and  oxidizer  yield  a  stoichiometric  mixture  fraction  Zst  =  0.143.  The  stoichiometric 
scalar  dissipation  rate  at  extinction  for  the  two  streams  (at  300  K  and  1  atm)  is  175  s-1, 
computed  in  a  counterflow  geometry  employing  the  reduced  mechanism.  The  background 
pressure,  initialized  at  1  atm,  increases  to  approximately  2  atm  towards  the  end  of  the  sim¬ 
ulation  due  to  heat  release.  Given  the  modest  increase  in  pressure  during  the  simulation, 
pressure  effects  are  not  expected  to  play  any  significant  role  in  the  physical  processes  of 
interest. 

The  reactive  scalar  fields  are  initialized  as  follows.  Temperature  and  chemical  species 
mass  fractions  are  taken  from  a  representative  one-dimensional  flamelet  solution  and  mapped 
from  mixture  fraction  space  onto  the  vertical  coordinate  according  to  a  specified  hyperbolic 
tangent  mixture  fraction  spatial  profile.  The  one-dimensional  flamelet  solution  is  obtained 
at  a  prescribed  stoichiometric  scalar  dissipation  rate  yst  =  60s-1,  corresponding  to  approxi¬ 
mately  one  third  of  the  extinction  conditions.  The  velocity  field  is  initialized  with  isotropic 
turbulence  fluctuations  of  a  prescribed  spectrum  and  Re>  =  170  at  the  start  of  the  simula¬ 
tion.  Additionally,  no  soot  is  present  in  the  domain  at  time  zero. 
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Figure  30:  Scalar  fields  of  temperature  (a,b)  and  scalar  dissipation  rate  in  logarithmic  scale 
(c,d)  at  1.5  and  10  ms  (top  to  bottom).  The  thick  solid  black  line  indicates  the  stoichiometric 
iso-contour,  while  the  thin  black  line  indicates  the  Z  =  0.3  iso-contour. 

Notwithstanding  the  two-dimensional  treatment  of  turbulence,  the  chosen  computational 
configuration  and  parameters  are  a  reasonable  approximation  of  the  dynamics  and  turbulence 
scales  characteristic  of  soot  formation  in  the  near  field  of  a  low  Reynolds  number  nonpremixed 
turbulent  jet  flame. 

4.2.2  Turbulent  Flame  Dynamics 

Due  to  the  strong  dependence  of  soot  growth  on  flame  parameters  such  as  temperature,  equiv¬ 
alence  ratio,  and  strain  rate,  turbulence/chemistry  interaction  is  investigated  first.  Turbulent 
fluid  motion  wrinkles  the  flame  and  causes  the  stoichiometric  contour  to  curve  and  stretch 
throughout  the  do-  main.  Instantaneous  snapshots  of  temperature  and  scalar  dissipation 
rate  at  times  1.5  and  10  ms  are  shown  in  Fig.  30. 

Scalar  dissipation  rate  is  distributed  spatially  in  thin,  elongated  layers.  A  few  locations 
along  the  flame  display  significant  weakening  and  extinction  due  to  high  scalar  dissipation 
rate  ( Xmax  ~  490s_1  at  1.5  ms).  For  instance,  local  extinction  can  be  observed  at  10  ms  in 
the  top  right  corner  of  the  domain,  where  the  stoichiometric  iso-contour  passes  through  a 
region  of  low  temperatures.  At  locations  where  flame  extinction  occurs,  the  mass  fraction  of 
OH  drops  from  its  peak  value  of  3.6  x  10-3  to  negligible  amounts  (not  shown). 
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Figure  31:  Scalar  fields  of  acetylene  (a,b)  and  naphthalene  (c,d)  at  1.5  and  10  ms  (top  to 
bottom).  The  thick  solid  black  line  indicates  the  stoichiometric  iso-contour,  while  the  thin 
black  line  indicates  the  Z  —  0.3  iso-contour. 
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In  the  soot  model  used  in  this  work,  two  species  have  significant  impact  on  soot  yield, 
namely  naphthalene  and  acetylene.  The  concentration  of  naphthalene  (CioHg)  governs  the 
rates  of  nucleation  and  condensation,  while  the  concentration  of  acetylene  (C2H2)  is  used  in 
evaluating  the  growth  rate  by  surface  reactions.  Figure  31  shows  acetylene  and  naphthalene 
mass  fractions  at  times  1.5  and  10  ms.  Compared  to  temperature,  the  mass  fractions  of 
naphthalene  and  acetylene  are  affected  more  significantly  by  turbulence  and  by  spatial  and 
temporal  variations  in  the  scalar  dissipation  rate.  At  10  ms,  the  mass  fraction  of  naphthalene 
is  highest  in  regions  of  low  scalar  dissipation  rate,  giving  rise  to  concentration  patches  away 
from  the  flame,  along  the  Z  =  0.3  iso-contour.  The  mass  fraction  of  acetylene  also  responds 
to  changing  values  of  scalar  dissipation  rate,  displaying  similar  spatial  patterns.  However, 
the  sensitivity  of  acetylene  to  turbulent  stretching  appears  to  be  less  pronounced  than  for 
naphthalene. 

To  better  analyze  the  effects  of  scalar  dissipation  rate  on  the  yield  of  acetylene  and 
naphthalene,  the  DNS  results  are  compared  to  the  solutions  of  the  steady  flamelet  equations 
and  are  shown  in  Fig.  32.  The  data  are  sampled  from  the  DNS  calculation  at  5  ms  and  along 
the  Z  =  0.3  iso-contour.  The  Z  =  0.3  location  was  chosen  as  it  corresponds  to  the  peak 
mass  fraction  for  both  acetylene  and  naphthalene  as  found  in  the  one- dimensional  flamelet 
calculations.  For  comparison,  the  one-dimensional  flamelet  calculations  were  performed  with 
the  same  assumptions  used  for  the  DNS  (e.g.  unity  Lewis  number,  no  Soret  and  Dufour 
effects,  etc.)  and  the  same  reduced  chemical  mechanism.  The  background  pressure  is  taken 
to  be  that  of  the  DNS  at  5  ms  (1.19  bar). 

Even  in  the  absence  of  turbulent  fluctuations,  the  mass  fraction  of  naphthalene  displays 
a  significantly  stronger  sensitivity  to  scalar  dissipation  rate  compared  to  acetylene.  With  the 
exception  of  the  data  points  at  low  values  of  scalar  dissipation  rates,  the  one-dimensional 
flamelet  solution  describes  well  the  dependence  of  the  mean  peak  mass  fraction  on  scalar 
dissipation  rate  shown  in  the  DNS  data. 

While  the  dependence  of  peak  mass  fraction  on  scalar  dissipation  rate  observed  in  the 
DNS  data  is  described  well  by  the  steady  flamelet  solution,  significant  scatter  exists  around 
the  mean  values.  The  magnitude  of  the  deviation  is  much  larger  for  naphthalene  than 
for  acetylene.  This  scatter  may  be  explained  by  the  rapidly  changing  turbulent  flow  field 
and  the  slowly  adjusting  chemical  species.  Furthermore,  the  flamelet  solution  at  low  scalar 
dissipation  rate  overestimates  the  mass  fraction  of  both  naphthalene  and  acetylene.  The 
difference  is  due  to  unsteady  flamelet  effects,  which  are  known  to  be  more  pronounced  at 
low  scalar  dissipation  rate  [83].  The  results  show  that  naphthalene  exhibits  much  stronger 
unsteady  effects  than  acetylene,  due  to  the  slower  chemistry  of  naphthalene  (lower  Damkohler 
number).  A  similar  behavior  is  expected  for  PAH  larger  than  naphthalene. 

4.2.3  Soot  Dynamics  and  Morphology 

The  emphasis  of  the  present  work  is  identifying  the  key  mechanisms  of  soot  growth  as  well 
as  characterizing  the  effects  of  turbulent  mixing  on  the  evolution  of  soot.  The  morphology 
of  soot  particles  is  investigated  also  by  focusing  on  aggregate  properties  such  as  the  primary 
particle  diameter  and  the  number  of  primary  particles. 
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X  at  Z  =  0.3  (1/s) 

Figure  32:  Mass  fraction  of  acetylene  and  naphthalene  as  a  function  of  scalar  dissipation 
rate.  Samples  taken  at  Z  —  0.3.  DNS  data  sampled  at  5  ms.  Acetylene  (top);  naphthalene 
(bottom):  —  steady  flamelet  solution;  •  ,  DNS  data. 
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4.2.3. 1  Soot  evolution  overview  Figure  33  shows  instantaneous  snapshots  of  soot 
number  density  and  soot  volume  fraction  fields  taken  at  10  and  15  ms.  The  spatial  distri¬ 
bution  of  soot  exhibits  two  distinct  patterns:  (a)  small  islands/patches  of  large  values  of 
soot  number  density  and  volume  fraction;  and  (b)  sharp  fronts  and  filament-like  patterns. 
While  the  islands  appear  to  accumulate  preferentially  along  the  Z  =  0.3  iso-contour,  the 
filaments  penetrate  deep  into  the  fuel  stream,  marking  the  edges  of  eddies.  The  number  and 
length  of  the  filaments  increase  with  time  as  turbulent  eddies  stretch  soot  patches  across 
the  domain  in  the  absence  of  diffusive  transport.  Due  to  the  high  Schmidt  number  (or  low 
mass  diffusivity)  typical  of  soot  transport,  soot  particles  do  not  diffuse,  and  their  transport 
is  controlled  mostly  by  convection  and  to  a  much  lesser  extent  by  thermophoresis. 

Soot  begins  nucleating  on  the  rich  side  of  the  flame  sheet  at  a  location  where  naphthalene 
abounds  (Z  ~  0.3).  From  the  onset  of  the  simulation  and  up  to  5  ms,  the  rate  of  increase 
of  soot  mass  is  dominated  by  the  nucleation  of  spherical  soot  particles.  At  early  times,  soot 
consists  of  a  large  number  of  small,  spherical  particles,  and  the  soot  particle  size  distribution 
is  monodisperse.  The  mean  diameter  and  number  of  primary  particles  at  1.5  ms  are  spatially 
homogeneous  and  equal  to  0.98  nm  and  1,  respectively  (not  shown).  Meanwhile,  naphthalene 
dimers  condense  on  the  surface  of  existing  particles.  During  the  simulation,  the  overall  rate 
of  condensation  increases  as  existing  particles  grow  and  more  surface  becomes  available. 
At  10  ms,  condensation  emerges  as  the  physical  process  contributing  most  to  soot  mass 
growth.  Intense  growth  occurs  in  regions  of  high  concentration  of  naphthalene  and  low 
scalar  dissipation  rate  as  is  evident  by  comparing  visually  Fig.  31(d)  and  Fig.  33(c).  For 
instance,  at  10  ms,  a  patch  of  naphthalene  mass  fraction  appears  near  the  bottom  left  corner 
of  Fig.  31(d)  where  there  is  a  peak  in  soot  volume  fraction  (see  Fig.  33(c)).  The  patch-like 
soot  patterns  observed  in  the  domain  are  identified  as  regions  of  naphthalene-based  soot 
growth  at  locations  of  low  scalar  dissipation  rate. 

The  patch-like  structure  of  soot  and  PAH  provides  insight  into  the  process  physics  and  its 
modeling.  From  a  physical  perspective,  the  sensitivity  of  PAH  species  to  hydrodynamic  strain 
and  unsteady  scalar  dissipation  rate  indicates  that  turbulent  mixing  in  technical  combustion 
devices  can  lead  to  soot  growth  suppression,  while  having  limited  effect  on  the  fuel  oxidation 
chemistry.  Based  on  the  results  presented,  the  dependence  of  soot  yield  on  the  mixing 
rate  can  be  explained  by  the  sensitivity  of  PAH  to  scalar  dissipation  rate.  Hydrodynamic 
strain  effects  are  expected  to  be  particularly  important  for  nonpremixed  and  low-temperature 
premixed  flames.  A  lesser  sensitivity  of  soot  growth  rates  to  strain  is  expected  for  high- 
temperature  premixed  flames  due  to  surface  reactions  being  the  dominant  growth  process. 
Surface  growth  rates  depend  on  the  concentration  of  acetylene,  which  is  less  sensitive  to 
strain.  From  a  modeling  perspective,  soot  growth  models  that  rely  on  smaller  hydrocarbon 
species  as  a  proxy  for  large  PAH  molecules  ignore  or  misrepresent  the  effects  of  turbulent 
mixing  and  hydrodynamic  strain  on  soot  formation  due  to  differences  between  the  Damkohler 
number  of  small  and  large  hydrocarbon  species.  For  example,  acetylene  is  often  used  as  an 
indicative  critical  species  in  the  soot  formation  process  in  order  to  simplify  the  treatment  of 
gas  phase  chemistry.  As  shown  in  this  work,  the  sensitivity  of  acetylene  to  turbulent  mixing 
is  significantly  less  pronounced  than  naphthalene.  Therefore,  the  application  of  a  soot  model 
relying  on  acetylene-based  mass  growth  to  a  configuration  for  which  PAH-based  growth  is 
the  dominant  growth  route  will  misrepresent  the  effect  of  hydrodynamic  strain  on  the  overall 
soot  yield. 
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Figure  33:  Scalar  fields  of  soot  particle  number  density,  N  (a,b)  and  soot  volume  fraction,  fv 
(c,d)  at  10  (top)  and  15  ms  (bottom).  The  thick  solid  black  line  indicates  the  stoichiometric 
iso-contour,  while  the  thin  black  line  indicates  the  Z  =  0.3  iso-contour. 
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Figure  34:  Conditional  statistics  of  the  displacement  speed  vz  along  the  Z  —  0.3  iso-contour 
versus  the  local  scalar  dissipation  rate  at  15  ms.  —  mean;  ----  standard  deviation;  •  data. 

4. 2. 3. 2  Soot  transport  and  growth  processes  Soot  growth  and  oxidation  rates  are  a 
function  of  temperature  and  mixture  composition.  While  local  mixture  fraction  determines 
nucleation,  con-  densation,  and  oxidation  rates,  it  provides  limited  information  about  quan¬ 
tities  such  as  soot  number  density  and  volume  fraction.  This  is  due  to  soot  experiencing 
significant  motion  in  mixture  fraction  space  because  of  differen-  tial  diffusion  effects  among 
soot  and  gas-phase  species.  Differential  diffusion  effects  are  important  in  soot  transport  and 
were  first  discussed  and  analyzed  by  Pitsch  et  al.  [84],  In  a  simple  candle  flame,  differ¬ 
ential  diffusion  causes  soot  to  be  transported  towards  the  flame  sheet  where  it  is  oxidized. 
This  process  is  similarly  responsible  for  soot  burnout  in  diesel  engines  and  aircraft  engine 
combustors. 

A  convenient  description  of  the  motion  of  a  mixture  fraction  (Z)  iso-surface  relative 
to  fluid  convection  is  provided  by  the  iso-surface  displacement  velocity  vz.  The  scalar 
quantity  vz  is  called  the  iso-surface  displacement  speed  and  takes  both  positive  and  negative 
values.  Due  to  negligible  mass  diffusivity,  soot  particles  are  convected  with  the  flow,  and 
the  displacement  speed  describes  the  movement  of  soot  relative  to  a  given  mixture  fraction 
iso-contour.  Thermophoresis  contributes  to  the  movement  of  soot  in  mixture  fraction  space 
also.  However,  for  the  conditions  of  the  present  DNS,  we  find  that  thermophoretic  velocities 
normal  to  the  Z  =  0.3  iso-contour  remain  well  below  3  cm/s  and  are  thus  neglected. 

Figure  34  shows  the  displacement  speed  vz  versus  scalar  dissipation  rate  xz  along  the 
Z  =  0.3  iso-contour,  which  corresponds  approximately  to  the  location  of  peak  soot  growth 
rates  and  soot  volume  fraction  for  the  chosen  stream  compositions.  The  data  are  sampled 
at  the  end  of  the  simulation  (15  ms).  Two  distinct  regimes  can  be  identified.  At  locations 
of  low  scalar  dissipation  rate,  the  absolute  value  of  the  displacement  speed  is  largest  and 
displays  both  negative  and  positive  values.  On  the  contrary,  in  regions  of  high  scalar  dissi¬ 
pation  rate,  the  displacement  speed  is  lowest  and  mostly  negative.  The  mean  displacement 
speed  conditioned  on  scalar  dissipation  rate  is  constant  across  the  range  of  scalar  dissipation 
rates  observed.  Mixing  reduces  the  peak  mixture  fraction  value  and  the  fraction  of  positive 
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Figure  35:  Statistics  of  the  curvature  and  normal  components  of  the  displacement  velocity. 
The  data  are  sampled  along  mixture  fraction  iso-contours  0.28  <  Z  <  0.32  at  15  ms.  The 
solid  lines  are  linear  least  square  fits  to  five  sets  of  data,  binned  according  to  the  value  of 
the  displacement  velocity:  (a)  160  <  vz  <  140,  (b)  60  <  vz  <  40,  (c)  10  <  vz  <  10, (d) 
40  <  vz  <  60,  (e)  80  <  vz  <  100  (cm/s).  For  each  set,  the  mean  values  of  the  two 
displacement  velocity  components  is  indicated  by  a  solid  square. 

displacement  speed  increases  as  soot  is  convected  towards  the  flame  sheet. 

In  order  to  better  understand  the  spatial  distribution  of  the  displacement  speed,  it  is 
written  as  the  sum  of  a  curvature  term  k,  proportional  to  the  iso-contour  curvature,  and  a 
normal  term,  proportional  to  the  spatial  derivative  of  the  mixture  fraction  gradient  along 
the  normal. 

Figure  35  illustrates  the  contribution  of  the  two  terms  contributing  to  the  overall  dis¬ 
placement  speed.  The  data  are  binned  with  respect  to  the  resulting  displacement  speed  to 
highlight  the  relative  importance  of  the  two  terms  for  a  given  value  of  vz  ■  Except  in  the 
case  of  very  small  values  of  the  displacement  speed,  the  curvature  and  normal  terms  have 
the  same  sign  and  act  synergistically.  This  fact  allows  to  judge  the  sign  of  the  displacement 
speed  from  the  sign  of  the  curvature  or  normal  component  alone.  Also,  the  mean  values  of 
the  two  components  are  similar,  indicating  that  (on  average)  they  contribute  equally  to  the 
resulting  displacement  speed.  Finally,  at  locations  where  the  curvature  term  is  the  largest 
in  the  domain,  the  normal  term  is  negligible.  Hence,  at  locations  of  maximum  iso-contour 
curvature,  the  displacement  speed  is  solely  due  to  the  curvature  component. 

Based  on  the  data  presented,  we  conclude  that  soot  particles  in  regions  of  positive  iso¬ 
contour  curvature  (center  of  curvature  on  the  lean  side  of  the  iso-contour)  move  towards  richer 
mixtures  (negative  displacement  speed);  the  opposite  occurs  for  soot  particles  in  regions  of 
negative  curvature  and  positive  displacement  speed. 

An  important  connection  between  soot  volume  fraction  and  iso-contour  curvature  is  es¬ 
tablished  by  the  correlation  between  curvature  and  scalar  dissipation  rate  along  the  Z  =  0.3 
iso-contour.  The  data  are  shown  in  Fig.  36  where  a  negative  correlation  between  curvature 
and  scalar  dissipation  rate  is  apparent:  regions  of  high  curvature  are  characterized  by  low 
scalar  dissipation  rate.  Due  to  the  sensitivity  of  PAH  to  scalar  dissipation  rate,  regions  of 
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Figure  36:  Statistics  of  scalar  dissipation  rate  conditioned  on  the  absolute  value  of  the  iso¬ 
contour  curvature  at  15  ms.  —  mean;  ----  standard  deviation;  •  data.  The  data  are  sampled 
along  the  Z  =  0.3  iso-contour. 

high  curvature  are  characterized  by  sustained  soot  nucleation  and  growth  and  (eventually) 
the  highest  soot  volume  fraction  (see  Fig.  37).  The  mean  soot  volume  fraction  along  three 
iso-contours  around  the  location  of  peak  growth  (Z  =  {0.23,  0.26,  0.3,  0.34})  increases  for  in¬ 
creasing  absolute  value  of  the  curvature,  reaching  a  minimum  at  zero  curvature.  The  mean 
soot  volume  fraction  profiles  are  approximately  symmetric  for  positive  and  negative  curva¬ 
tures.  A  small,  but  obvious  bias  in  favor  of  positive  curvature  and  negative  displacement 
speed  appears  as  the  data  are  sampled  along  increasingly  rich  iso-contours. 

Due  to  the  symmetry  in  the  statistics  of  the  displacement  speed  at  low  scalar  dissipa¬ 
tion  rate,  movement  towards  leaner  or  richer  compositions  is  almost  equally  likely.  The 
displacement  direction  established  by  the  sign  of  the  iso-contour  curvature  determines  the 
subsequent  soot  evolution.  Soot  moving  towards  the  flame  (negative  curvature)  first  grows 
by  surface  reaction  via  the  HACA  mechanism  at  locations  of  high  temperature  and  acetylene 
concentration.  As  the  displacement  continues,  soot  is  oxidized  near  the  flame  sheet.  The 
occurrence  of  surface  growth  followed  by  oxidation  at  the  flame  sheet  was  first  described  by 
Pitsch  et  al.  [84]  in  the  context  of  unsteady  flamelets  and  recently  demonstrated  in  two  DNS 
studies  of  turbulent  sooting  flames  [28,  29]. 

Due  to  the  inverse  correlation  between  curvature  and  scalar  dissipation  rate,  the  most 
rapid  displacement  in  mixture  fraction  space  occurs  at  locations  where  abundant  PAH  species 
generate  the  highest  values  of  soot  volume  fraction.  The  correlation  between  high  values 
of  soot  growth  rates  and  pronounced  differential  diffusion  motion  in  Z-space  explains  the 
observed  soot  distribution  and  growth  patterns  in  mixture  fraction  space.  First,  soot  nu¬ 
cleates  around  Z  =  0.26  in  regions  of  low  scalar  dissipation  rate.  Soot  particles  in  those 
regions  characterized  by  positive  curvature  are  displaced  towards  richer  mixtures.  During 
this  relative  motion,  mass  growth  continues  sustained  by  PAH  species  condensation  on  soot 
particles.  Growth  by  condensation  eventually  surpasses  growth  by  nucleation  as  more  soot 
surface  becomes  available.  Since  the  concentration  of  naphthalene  decreases  monotonically 
for  Z  >  0.26,  soot  growth  rates  decrease  as  the  aggregates  are  displaced  further  towards  richer 
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Figure  37:  Mean  soot  volume  fraction  along  four  mixture  fraction  iso-contours  (Z  = 
{0.23,  0.26,  0.3,  0.34})  conditioned  on  the  local  value  of  the  iso-contour  curvature.  Data 
sampled  at  15  ms. 

mixtures  and  cross  the  Z  =  0.3  iso-contour.  Soot  volume  fraction  peaks  approximately  at 
Z  =  0.31,  i.e.  at  richer  conditions  than  those  of  peak  nucleation  and  peak  condensation  due 
to  a  drift  in  mixture  fraction  space  caused  by  differential  diffusion.  On  average,  soot  growth 
due  to  surface  reactions  is  maximum  at  Z  =  0.23,  i.e.  at  mixture  compositions  leaner  than 
those  characteristic  of  peak  nucleation. 

Despite  the  movement  of  soot  in  mixture  fraction  space  being  nearly  symmetric,  the 
spatial  distribution  of  soot  volume  fraction  is  significantly  biased  towards  the  fuel  stream 
due  to  higher  PAH-based  growth  rates  at  rich  mixtures  and  an  oxidizing  environment  next 
to  the  flame  sheet. 

4.3  LES  Model  Validation  and  Simulation  Results 
4.3.1  Validation  Case:  Delft  Flame  III 

The  configuration  that  will  be  used  to  validate  the  model  is  Delft  Flame  III  [85],  a  piloted 
turbulent  nonpremixed  jet  flame  fueled  by  commercial  natural  gas.  At  the  burner  exit, 
a  6  mm  diameter  fuel  jet  is  surrounded  by  a  primary  annular  coflow  of  air  with  an  inner 
diameter  of  15  mm,  which  decreases  from  30  mm  over  the  last  60  mm  of  the  burner,  and  an 
outer  diameter  of  45  mm.  The  flame  is  stabilized  by  twelve  small  pilot  flames  arranged  in 
a  ring  which  emanate  from  the  rim  separating  the  fuel  stream  from  the  primary  air.  These 
pilot  flames  are  part  of  an  insert  in  the  fuel  tube  which  causes  the  fuel  jet  diameter  to 
decrease  from  8  mm  at  a  location  16  mm  upstream  of  the  burner  exit.  The  bulk  velocity  of 
the  fuel  stream  is  21.9  m/s  giving  a  jet  Reynolds  number  of  roughly  8200  based  on  the  fuel 
jet  diameter  and  bulk  velocity.  The  bulk  velocity  of  the  primary  air  stream  is  4.4  m/s,  and 
the  coflow  velocity  is  0.3  m/s. 

Several  experimental  campaigns  have  been  carried  out  with  the  burner,  and  three  of  these 
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datasets  will  be  used  for  validation  in  this  work.  Two-component  velocity  measurements  were 
made  with  Laser  Doppler  Anemometry  (LDA)  by  Stroomer  [1],  Subsequent  measurements 
of  the  thermochemical  scalars  were  made  by  Nooren  et  al.  [2]  using  the  Raman-Rayleigh-LIF 
(Laser  Induced  Fluorescence)  technique.  These  measurement  were  confined  to  a  region  closer 
to  the  burner  where  soot  is  not  present.  Further  downstream,  the  soot  volume  fraction  was 
measured  using  Laser  Induced  Incandescence  (LII)  by  Qamar  et  al.  [3]. 

Since  the  measurements  were  made  in  three  different  places,  the  composition  of  com¬ 
mercial  natural  gas  was  not  exactly  the  same.  In  the  latter  two  studies  [2,  3],  the  natural 
gas  was  diluted  with  N2  to  match  the  adiabatic  flame  temperature  of  the  first  study.  In 
this  work,  the  fuel  composition  from  the  soot  volume  fraction  measurements  will  be  used. 
However,  analysis  has  shown  that  the  results  of  this  work  are  not  sensitive  to  a  specific  fuel 
composition. 

4.3.2  Computational  Implementation 

The  models  presented  in  the  previous  section  are  implemented  in  NGA,  a  structured  finite 
difference  code  for  low  Mach  number  turbulent  reacting  flows  that  was  initially  developed  at 
Stanford.  The  code  is  based  on  the  numerical  methods  of  Desjardins  et  al.  [86].  In  this  work, 
the  momentum  and  continuity  equations  are  discretized  with  second  order  spatial  operators. 
The  scalar  equations  required  in  the  presumed  PDF  approach  are  spatially  discretized  with 
a  bounded  QUICK  scheme  [87].  All  of  the  subfilter  stresses  and  fluxes  are  closed  with  the 
dynamic  model  [88-90]  with  Lagrangian  averaging  [91,  92],  In  order  to  accommodate  the 
large  four- dimensional  tables  required  for  the  RFPV  and  soot  models,  the  code  has  been 
parallelized  using  a  hybrid  method  with  one  Message  Passing  Interface  (MPI)  process  per 
cluster  node  and  several  OpenMP  threads  per  node. 

In  order  to  generate  the  boundary  conditions  for  the  simulation,  two  additional  smaller 
simulations  are  first  performed.  First,  a  fully  developed  pipe/ annulus  flow  is  computed  for  the 
burner  upstream  of  the  exit  region.  Second,  the  results  of  this  simulation  are  used  as  inflow 
conditions  for  a  simulation  of  the  last  few  centimeters  of  the  burner  where  the  geometries 
of  the  primary  air  annulus  and  fuel  jet  change.  The  results  of  this  second  simulation  are 
then  used  as  inflow  conditions  for  the  main  flame  simulation.  The  secondary  air  co-flow  is 
prescribed  as  a  bulk  flow.  The  pilot  geometry  is  modeled  as  an  annulus  rather  than  twelve 
individual  streams  to  preserve  axisymmetry,  and  the  inflow  velocity  is  specified  as  a  bulk 
flow.  The  pilot  velocity  is  reduced  (since  the  pilot  area  is  larger)  such  that  the  mass  flow  rate 
matches  the  experiment  and  was  found  to  be  sufficiently  strong  to  keep  the  flame  attached  to 
the  burner.  The  reduced  pilot  velocity  has  an  additional  benefit  of  allowing  a  large  timestep 
to  be  taken,  which  is  absolutely  vital  since  soot  in  the  downstream  region  of  the  flame  evolves 
very  slowly. 

The  computational  domain  is  discretized  with  384,  192,  and  64  points  in  the  axial,  radial, 
and  circumferential  directions,  respectively.  The  domain  is  decomposed  with  96  processors 
with  a  total  computational  cost  of  about  8.5  s//rs.  Statistics  are  collected  over  a  period  of 
approximately  O  (450  ms)  for  a  total  computational  cost  of  roughly  O  (100,000  cpu-hrs). 
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Figure  38:  Radial  profiles  of  the  Favre-averaged  axial  velocity.  The  axial  locations  are  indi¬ 
cated.  Lines  are  the  LES,  and  symbols  are  the  experimental  measurements  of  Stroomer  [1], 


Figure  39:  Radial  profiles  of  the  Favre-averaged  mixture  fraction.  The  axial  locations  are 
indicated.  Lines  are  the  LES,  and  symbols  are  the  experimental  measurements  of  Nooren  et 
al  [2], 


4.3.3  Presumed  PDF  Method  Results 

4.3.3. 1  Upstream  Results:  Velocity  and  Scalars  In  this  section,  the  LES  results  are 
compared  to  the  experimental  velocity  and  scalar  measurements  upstream  where  soot  is  not 
present  [1,  2].  In  addition  to  the  soot  model,  this  work  is  also  the  first  LES  of  Delft  Flame 
III.  Where  appropriate,  results  of  previous  studies  of  Delft  Flame  III  utilizing  RANS  will 
be  discussed.  While  accurate  predictions  of  the  velocity  and  scalars  upstream  do  not  imply 
that  these  quantities  are  predicted  accurately  downstream  where  soot  is  present,  accurate 
predictions  upstream  do  increase  the  confidence  in  the  corresponding  downstream  quantities 
which  cannot  be  validated  directly  against  experimental  measurements. 

Radial  profiles  of  the  axial  velocity  are  shown  in  Fig.  38.  Overall,  agreement  between  the 
LES  and  the  experimental  measurements  is  excellent.  At  the  furthest  station  downstream, 
the  LES  predicts  a  slightly  smaller  velocity  than  the  experimental  measurements,  but  this 
small  deviation  is  certainly  within  the  experimental  uncertainty  [1],  In  a  previous  RANS 
study  [93] ,  the  axial  velocity  at  the  furthest  downstream  station  was  significantly  underpre¬ 
dicted  for  nearly  all  of  the  turbulence  and  combustion  models  investigated.  Compared  to 
this  previous  study,  the  present  LES  predicts  the  axial  velocity  much  more  accurately. 

Radial  profiles  of  the  mixture  fraction  are  shown  in  Fig.  39.  Overall,  like  the  velocity,  the 
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Figure  40:  Radial  profiles  of  the  Favre-averaged  temperature.  The  axial  locations  are  in¬ 
dicated.  Lines  are  the  LES,  and  symbols  are  the  experimental  measurements  of  Nooren  et 
al  [2], 


LES  agrees  very  well  with  the  experimental  measurements.  The  total  potential  experimental 
uncertainty  for  the  mixture  fraction  is  just  less  than  10%  (though  possibly  up  to  nearly  25% 
in  rich  mixtures)  [2],  and  most  of  the  deviations  from  the  experimental  measurements  are 
within  this  uncertainty.  However,  there  are  two  minor  discrepancies  worth  noting.  First, 
the  spreading  rate  of  the  jet  at  the  first  station  near  the  nozzle  is  slightly  too  fast.  Second, 
at  the  two  downstream  stations,  the  centerline  mixture  fraction  is  slightly  overpredicted. 
These  minor  differences  become  more  apparent  in  the  temperature  and  species  mass  fractions 
profiles  that  follow. 

In  previous  RANS  studies  [93-96],  the  same  two  discrepancies  are  observed,  that  is, 
an  initial  jet  spreading  rate  that  is  too  fast  and  an  overprediction  of  the  centerline  mixture 
fraction  at  downstream  locations.  However,  in  the  previous  RANS  studies,  an  overprediction 
of  the  jet  spreading  rate  remains  even  at  the  stations  furthest  downstream  from  the  nozzle. 
In  the  LES,  the  spreading  rate  of  the  jet  seems  to  be  only  slightly  overpredicted  very  close 
to  the  burner. 

The  minor  discrepancy  in  the  jet  spreading  rate  near  the  burner  is  a  consequence  of 
the  model  for  the  subfilter  mixture  fraction  variance.  As  discussed  previously,  an  algebraic 
model  cannot  be  used  in  the  downstream  region  of  the  flame  where  soot  is  present.  However, 
upstream,  the  dynamic  algebraic  model  of  Pierce  and  Moin  [54]  can  be  used.  Using  this 
model  (results  not  shown),  the  mixture  fraction  variance  was  significantly  smaller,  consistent 
with  the  work  of  Kaul  and  coworkers  [55,  56],  and  the  initial  jet  spreading  rate  was  well 
captured.  Presumably,  the  correct  spreading  rate  could  also  be  captured  with  the  transport 
equation  model  with  the  use  of  a  more  accurate  model  for  the  subfilter  scalar  dissipation 
rate.  However,  it  is  worth  noting  that  the  initial  spreading  rate  was  not  sensitive  to  the 
value  of  the  coefficient  in  Eq.  31,  so  a  change  in  the  model  form  would  be  needed. 

Radial  profiles  of  the  temperature  are  shown  in  Fig.  40.  Overall,  both  the  shapes  and 
magnitudes  of  the  experimental  profiles  are  predicted  relatively  accurately.  Close  to  the 
burner,  the  slightly  broadened  temperature  profile  is  consistent  with  the  the  corresponding 
mixture  fraction  profile  in  Fig.  39.  Furthest  from  the  burner,  the  temperature  is  somewhat 
overpredicted  near  the  centerline.  However,  the  experimental  temperatures  obtained  with 
the  Raman-Rayleigh-LIF  technique  may  underpredict  the  temperature  by  up  to  70  K,  es- 
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Species 

Uncertainty  (%) 
Lean  Rich 

CO 

20 

25 

co2 

12 

8 

H2 

17 

37 

h2o 

11 

18 

Table  1:  Total  potential  experimental  uncertainty  for  the  species  mass  fraction  measure¬ 
ments  [2], 

pecially  closer  to  the  centerline  where  CH4  interferes  with  the  temperature  measurement. 
Subsequent  temperature  measurements  using  CARS  [97]  confirmed  the  underprediction  of 
the  temperature  with  the  Raman-Rayleigh-LIF  technique. 

In  the  previous  RANS  studies  utilizing  an  adiabatic  steady  flamelet  model  [93],  the  tem¬ 
perature  was  globally  overpredicted.  Habibi  et  al.  [96]  considered  radiation  effects  with  a 
steady  flamelet  model.  However,  while  the  inclusion  of  radiation  improved  the  prediction 
of  the  temperature  at  the  measurement  station  furthest  from  the  burner,  the  temperature 
closer  to  the  burner  was  overpredicted.  The  overprediction  is  a  consequence  of  using  only 
burning  flamelets  with  the  steady  flamelet  model  which  does  not  capture  the  minor  local 
extinction  in  the  flame.  With  transported  PDF  methods  [93-95],  local  extinction  can  be 
captured  with  RANS  models.  However,  the  predicted  temperature  is  very  sensitive  to  the 
amount  of  local  extinction  predicted  by  the  mixing  model  and  the  value  of  any  specified 
coefficient.  The  present  LES  RFPV  model  is  an  improvement  upon  both  classes  of  RANS 
combustion  models,  capturing  the  correct  amount  of  local  extinction  and  accurately  predict¬ 
ing  the  temperature  without  a  strong  sensitivity  to  a  constant  that  must  be  specified  for  a 
mixing  model.  Admittedly,  a  constant  must  be  specified  for  the  subfilter  scalar  dissipation 
rate  model  in  Eq.  31,  but  the  predicted  temperatures  at  the  measurement  stations  considered 
in  Fig.  40  are  not  sensitive  to  this  constant. 

Radial  profiles  of  the  major  products  of  combustion  are  shown  in  Fig.  41.  Nearly  all 
of  the  discrepancies  observed  for  the  species  are  consistent  with  discrepancies  observed  for 
other  quantities.  For  example,  CO  and  H2  are  overpredicted  near  the  centerline  where  the 
mixture  fraction  is  overpredicted.  Considering  the  larger  experimental  uncertainties  for  the 
individual  species  mass  fractions,  which  are  summarized  in  Table  1,  overall  agreement  with 
the  experimental  measurements  can  be  considered  good. 

As  discussed  above,  near  the  burner  where  soot  is  not  present,  the  subfilter  variance 
could  also  be  modeled  with  an  algebraic  model.  With  an  algebraic  model,  both  CO  and 
H2  show  marked  decreases  (not  shown)  compared  to  the  transport  equation  model,  bringing 
them  into  better  agreement  with  the  experimental  measurements.  In  addition,  just  as  form 
the  mixture  fraction,  the  species  do  not  show  a  significant  sensitivity  to  the  coefficient  in 
the  subfilter  scalar  dissipation  rate  model,  at  least  compared  to  the  sensitivity  of  the  species 
to  the  subfilter  variance  model.  This  observation  further  highlights  the  need  for  improved 
subfilter  scalar  dissipation  rate  models. 

4. 3. 3. 2  Downstream  Results:  Soot  In  this  section,  the  downstream  LES  results  are 
presented.  First,  the  evolution  of  soot  in  the  jet  flame  is  discussed  qualitatively.  Particular 
attention  will  be  paid  to  the  relative  roles  of  the  various  formation,  growth,  and  destruction 
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Figure  41:  Radial  profiles  of  the  Favre-averaged  mass  fractions.  The  axial  locations  are 
indicated.  Lines  are  the  LES,  and  symbols  are  the  experimental  measurements  of  Nooren  et 
al  [2], 
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Figure  42:  Instantaneous  soot  properties.  The  solid  black  lines  correspond  to  the  stoichio¬ 
metric  mixture  fraction  iso-contour.  The  horizontal  lines  demarcate  intervals  of  25x/D  from 


25  to  125. 


processes.  Then,  the  LES  results  are  compared  to  the  experimental  soot  volume  fraction 
measurements  [3].  Finally,  soot  is  sensitive  to  the  coefficient  in  subfilter  scalar  dissipation 
rate  model  (Eq.  31),  and  this  sensitivity  is  investigated  in  detail. 

The  instantaneous  fields  of  a  few  soot  properties  are  shown  in  Fig.  42.  Several  qualitative 
features  are  immediately  apparent.  First,  soot  is  confined  only  to  regions  of  rich  mixture 
fraction.  When  soot  reaches  stoichiometric  mixtures,  oxidation  by  OH  proceeds  so  quickly 
that  soot  never  has  an  appreciable  residence  time  at  lean  mixtures.  Second,  soot  is  not 
homogeneous  at  a  given  mixture  fraction  or,  in  other  words,  is  spatially  intermittent,  which 
is  a  consequence  of  the  sensitivity  of  PAH  to  the  local  scalar  dissipation  rate.  Spatial 
intermittency  also  arises  from  the  fact  that  soot  has  no  molecular  diffusivity.  However,  most 
of  this  structure  is  lost  in  LES  due  to  the  implicit  spatial  filtering.  Conversely,  temporal 
intermittency  is  readily  observed  with  more  than  a  factor  of  three  variation  in  the  maximum 
soot  volume  fraction  in  the  domain.  Finally,  the  largest  soot  volume  fractions  tend  to 
be  present  at  mixture  fractions  significantly  greater  than  stoichiometric.  This  fact  is  a 
consequence  of  the  active  soot  growth  mechanisms  and  will  be  discussed  in  greater  detail  in 
the  subsequent  paragraphs. 

Each  of  the  four  regions  indicated  in  Fig.  42  corresponds  to  a  distinct  phase  of  the  soot 
evolution.  In  the  first  region  (x/D  <  50),  the  soot  number  density  is  large  while  the  volume 
fraction  and  average  primary  particle  diameter  are  relatively  small.  Here,  nucleation  of  soot 
particles  from  PAH  is  the  lone  active  process,  resulting  in  many  small  particles.  At  slightly 
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Figure  43:  Time- averaged 


properties.  All  of  the 


(c)  Primary  Particle  Diameter, 
dp  [nm] 


are  the  same  as  in  Fig.  42. 


larger  heights  above  the  burner  (50  <  x/D  <  75),  the  number  density  remains  roughly 
constant  while  the  volume  fraction  and  average  primary  particle  diameter  begin  to  increase. 
The  particles  formed  in  the  lower  region  are  now  numerous  enough  that  some  of  the  PAH 
begins  to  condense  onto  existing  particles  rather  than  forming  new  particles.  At  higher 
heights  yet  (75  <  x/D  <  100),  the  volume  fraction  and  primary  particle  diameter  attain 
their  maximum  values,  and  the  number  density  begins  to  decrease.  At  these  locations,  nearly 
all  of  the  PAH  condenses  onto  the  now  larger  particles,  and  coagulation  results  in  fewer  but 
larger  aggregates.  Finally,  far  downstream  (x/D  >  100),  regions  of  rich  mixture  fraction 
start  to  disappear  as  the  flame  closes,  and  the  particles  are  oxidized  by  OH  as  they  approach 
and  finally  pass  through  the  stoichiometric  mixture  fraction  iso-contour. 

Just  as  observed  in  the  DNS  study  of  Bisetti  et  al.  [24],  surface  growth  is  found  to  play 
a  limited  role  in  the  evolution  of  soot.  In  this  nonpremixed  flame,  soot  nucleates  around  a 
mixture  fraction  of  0.2.  Then,  according  to  the  DNS  study,  transport  in  mixture  fraction 
space  toward  leaner  or  richer  mixtures  is  roughly  equally  probable.  If  the  soot  is  transported 
to  richer  mixture  fraction,  it  lingers  and  grows  according  to  the  processes  outlined  above 
(condensation  followed  by  coagulation).  If  the  soot  is  transported  to  leaner  mixture  fractions, 
it  first  grows  by  surface  growth  but  is  then  immediately  oxidized  as  it  passes  through  the 
stoichiometric  mixture  fraction  iso-contour.  However,  since  these  particles  are  small,  the 
cumulative  growth  by  surface  growth  is  relatively  small,  and,  as  a  result,  surface  growth 
does  not  contribute  significantly  to  the  total  soot  formation  rate  in  this  flame. 

The  time-averaged  fields  of  a  few  soot  properties  are  shown  in  Fig.  43.  The  variation  in  the 
quantities  moving  downstream  from  the  burner  is  very  similar  to  the  instantaneous  images 
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Figure  44:  Axial  profiles  of  the  averaged  soot  volume  fraction  and  total  intermittency.  Lines 
are  the  LES,  and  symbols  are  the  experimental  measurements  of  Qamar  et  al.  [3].  The  solid 
red  line  is  the  nominal  constant  for  the  subfilter  dissipation  rate  model  (Eq.  31);  the  dashed 
green  line  is  half  the  constant;  and  the  dotted  blue  line  is  double  the  constant.  Note  that 
for  the  total  intermittency  the  LES  results  and  experimental  measurements  are  plotted  with 
different  vertical  scales. 


in  Fig.  42.  However,  note  that  the  maximum  time-averaged  volume  fraction  and  number 
density  are  significantly  lower  than  the  instantaneous  maxima.  This  is  a  direct  consequence 
of  the  intermittency  of  soot  as  a  result  of  the  sensitivity  of  PAH  to  the  scalar  dissipation 
rate  and  a  lesser  extent  to  the  zero  diffusivity  of  soot  since  the  subfilter  diffusivity  dominates 
in  LES.  Also,  in  Fig.  43,  significant  soot  volume  fraction  is  present  at  lean  mixtures  due  to 
the  turbulent  fluctuations  in  the  radial  direction  (and  axial  direction  at  the  flame  tip)  of 
the  stoichiometric  mixture  fraction  iso-contour.  With  a  RANS  model,  this  effect  would  not 
be  observed,  and  soot  would  be  confined  entirely  to  regions  of  rich  mixture  fraction  since  it 
would  be  oxidized  near  the  stoichiometric  mixture  fraction  iso-contour. 

The  axial  profile  of  the  mean  soot  volume  fraction  is  shown  in  Fig.  44  and  compared 
with  the  experimental  measurements  of  Qamar  et  al.  [3].  The  LES  model  is  shown  to 
reasonably  predict  the  magnitude  of  the  soot  volume  fraction,  although  this  measure  is 
sensitive  to  the  subfilter  scalar  dissipation  rate  model  as  discussed  in  the  subsequent  section. 
However,  soot  is  formed  much  earlier  in  the  simulation  than  in  the  experiment.  Subsequent 
analysis  revealed  the  onset  of  soot  was  not  sensitive  the  subfilter  dissipation  rate  model. 
In  this  region  only  nucleation  is  active,  which  depends  essentially  only  on  the  chemistry. 
Therefore,  without  additional  validation  data,  this  discrepancy  must  be  attributed  largely 
to  the  significant  uncertainty  in  the  kinetic  mechanism  for  PAH  formation  and  growth, 
particularly  for  methane  flames. 

The  axial  profile  of  the  total  intermittency  is  also  shown  in  Fig.  44  and  compared  with 
the  experimental  measurements.  The  total  intermittency  is  defined  by  Qamar  et  al.  [3]  as  the 
probability  of  the  soot  volume  fraction  being  less  than  0.1  ppb  (the  experimental  threshold) 
at  a  given  location  at  a  given  time.  In  the  LES,  this  measure  is  a  function  of  not  only 
the  volume  fraction  but  also  the  subfilter  intermittency.  While  the  LES  predicts  a  much 
lower  intermittency  than  the  experiments,  the  comparison  is  difficult  given  the  discrepancy 
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Figure  45:  Radial  profiles  of  soot  volume  fraction,  temperature,  and  mixture  fraction  in  the 
region  of  peak  soot  volume  fraction.  The  different  lines  are  the  same  as  in  Fig.  44. 


in  the  volume  fraction.  However,  like  the  experiments,  the  LES  shows  an  increase  in  the 
intermittency  in  the  oxidation  region  of  the  flame.  As  discussed  by  Qamar  et  al.  [3],  this 
indicates  a  reduction  in  frequency  of  soot  events  and  not  the  intensity  of  these  events  in  this 
region  of  the  flame.  Based  on  the  LES  results,  this  is  a  direct  consequence  of  rich  pockets 
of  mixture  fraction  intermittently  pinching  off  from  the  continuous  region  of  rich  mixture 
fraction  upstream.  This  phenomenon  can  be  observed  in  Fig.  42. 

As  alluded  to  previously,  the  LES  results  are  quites  sensitive  to  the  subfilter  dissipation 
rate  model.  Fig.  44  includes  results  for  not  only  the  nominal  value  of  the  constant  in 
Eq.  31  but  also  the  results  with  the  coefficient  doubled  and  halved.  As  shown  in  Fig.  44, 
the  maximum  time-averaged  volume  fraction  decreases  by  nearly  a  factor  of  two  when  the 
constant  in  Fig.  31  is  halved.  On  the  other  hand,  the  maximum  volume  fraction  only  increases 
by  about  25%  when  the  constant  is  doubled.  The  value  of  the  constant  does  not  appreciably 
affect  the  location  of  the  onset  of  soot,  the  flame  length,  or  the  intermittency. 

In  order  to  determine  the  reason  for  the  sensitivity  the  subfilter  dissipation  rate,  radial 
profiles  of  the  volume  fraction,  temperature,  and  mixture  fraction  are  shown  in  Fig.  45  at 
two  axial  locations  where  the  soot  volume  fraction  is  large.  At  both  locations,  there  is  only 
a  minor  difference  between  the  temperature  and  mixture  fraction  profiles,  some  of  which  is 
due  to  a  lack  of  complete  statistical  convergence,  but  the  soot  volume  fraction  profiles  are 
markedly  different.  This  indicates  that  the  subfilter  dissipation  rate  model  does  not  affect 
soot  through  the  background  flow  field  and  thermochemical  state  but  instead  through  the 
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Figure  46:  Conditional  soot  source  terms.  Lines  are  the  same  as  in  Fig.  44.  Positive 
values  correspond  to  the  combined  PAH  formation  and  growth  pathways  (nucleation  and 
condensation),  and  negative  values  correspond  to  oxidation. 


mixture  fraction  variance  itself. 

When  the  subfilter  dissipation  is  decreased  and  the  subfilter  variance  increased,  the  sub¬ 
filter  distribution  of  mixture  fraction  includes  both  leaner  values  and  richer  values.  Increased 
amounts  of  lean  mixture  fractions  increase  the  rate  of  oxidation,  and  as  a  result,  as  the  vari¬ 
ance  increases,  the  rate  of  oxidation  increases.  Conversely,  increased  amounts  of  rich  mixture 
fractions  increase  the  combined  rates  of  nucleation  and  condensation.  However,  the  effect  for 
rich  mixture  fractions  is  diminished  as  the  PAH  concentration  is  constrained  to  a  relatively 
narrow  profile  centered  around  approximately  Z  =  0.2.  This  effect  is  readily  seen  in  plots 
of  the  conditional  source  terms  in  Fig.  46.  For  larger  variances,  oxidation  occurs  at  larger 
mixture  fractions.  On  the  other  hand,  nucleation  and  condensation  are  virtually  unaffected. 
In  the  end,  for  the  smaller  subfilter  dissipation,  formation  and  growth  are  more  co-located 
with  oxidation  in  mixture  fraction  space,  and  the  total  soot  yield  is  reduced. 

An  additional  interesting  trend  in  Fig.  44  is  that  the  beginning  of  the  increase  in  the  soot 
volume  fraction  is  largely  independent  of  the  model  constant  for  the  subfilter  dissipation 
rate.  In  this  region,  the  centerline  line  mixture  fraction  is  still  sufficiently  large  that,  even 
for  large  subfilter  variances,  virtually  no  lean  mixtures  are  present  in  the  subfilter  distribu¬ 
tion  of  the  mixture  fraction.  As  a  result,  the  argument  above  for  increased  oxidation  with 
increased  subfilter  variance  does  not  apply.  Therefore,  the  premature  increase  in  the  volume 
fraction  compared  to  the  experiments  is  likely  due  to  the  significant  uncertainty  in  the  ki¬ 
netic  mechanism  for  PAH  formation  and  growth.  A  priori  analysis  of  the  flamelet  solutions 
revealed  significant  sensitivity  to  certain  reaction  pathways  (such  as  that  discussed  in  the 
previous  section),  which  could  easily  account  for  the  observed  discrepancies. 

4.3.4  Transported  PDF  Method  Results 

4.3.4. 1  Instantaneous  fields  Figure  47  shows  the  instantaneous  fields  of  temperature, 
PAH  mass  fraction,  soot  volume  fraction,  and  soot  number  density.  As  expected,  soot  forma¬ 
tion  is  predominant  in  the  fuel-rich  region  inside  the  stoichimetric  surface  where  moderately 
high  temperatures  aid  the  formation  of  PAH  molecules.  While  PAH  formation  is  initiated 
far  upstream  (x/d  >  20),  significant  soot  volume  fraction  is  observed  only  after  x/d=60. 
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Figure  47:  Instantaneous  contours  of  temperature,  PAH  mass  fraction,  soot  volume  fraction, 
and  soot  number  density.  Iso- value  of  stoichiometric  mixture  fraction  is  indicated  by  a  solid 
black  line. 


The  lag  between  PAH  growth  and  soot  formation  shows  the  time  scale  differences  between 
the  two  processes.  Similar  to  the  experiment,  soot  formation  is  found  to  be  highly  intermit¬ 
tent  with  sporadic  bursts  of  high  soot  volume  fraction  regions  followed  by  extended  periods 
of  low  soot  volume  fraction.  Prior  studies  have  demonstrated  that  soot  precursor  growth, 
especially  the  evolution  of  naphthalene,  is  highly  sensitive  to  the  strain  rates  in  the  flow. 
Consequently,  the  turbulent  features  that  a  fluid  packet  experiences  strongly  dictates  soot 
nucleation  [24],  Figure  47  also  shows  that  the  soot  number  density  starts  to  increase  far 
upstream  due  to  nucleations  up  to  x/d=20.  Further  downstream,  the  number  density  does 
not  decrease  substantially,  indicating  that  soot  coagulation  is  almost  balanced  by  persisting 
nucleation  process. 

4. 3. 4. 2  Soot  and  gas-phase  statistics  Figure  48  shows  the  time-averaged  mixture 
fraction  and  temperature  statistics  for  the  flame,  measured  far  upstream  of  the  region  of 
high  soot  volume  fraction.  Results  indicate  that  the  LES/PDF  predicts  the  flame  statistics 
reasonably  accurately,  including  the  RMS  statistics  of  the  flow.  In  the  PDF  technique, 
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x/d=41 .7 


x/d=41 .7 


Figure  48:  Comparison  of  time  averaged  mean  and  RMS  of  mixture  fraction  and  temper¬ 
ature  with  the  experiment  of  Nooren,  et.  al.  [2],  Simulation:  mean  ),  RMS  (----  ); 
Experiment:  mean  (■  ),  RMS  (A). 


Figure  49:  Comparison  of  soot  volume  fraction  from  the  simulation  with  the  experiment  of 
Qamar,  et.  al.  [3].  The  data  is  plotted  along  the  centerline. 


the  subfilter  mixing  model  (Eq.  45)  is  a  source  of  spatial  numerical  diffusion,  since  it  mixes 
particles  that  are  spatially  distributed  within  a  single  control  volume.  As  the  dissipation  rate 
increases,  the  higher  mixing  rates  lead  to  increased  numerical  diffusion  that  will  manifest  as 
reduced  temporal  fluctuations.  The  accuracy  of  the  RMS  statistics  provides  confidence  that 
the  PDF  numerics  have  not  unduly  affected  the  flame  computation. 

Figure  49  shows  the  soot  volume  fraction  obtained  from  the  LES/PDF  approach  and 
the  experiments.  The  simulation  results  over  predict  the  soot  volume  fraction  substantially, 
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and  also  show  an  earlier  inception  and  growth  compared  to  the  experiments.  Mueller  and 
Pitsch  [98]  reported  that  the  soot  profiles  showed  very  large  variations  with  the  dissipation 
rate  models  used.  Here,  a  dynamic  procedure  is  used  to  evaluate  the  dissipation  rate  [66]. 
While  this  approach  was  shown  to  correctly  predict  extinction  levels  in  partially-premixed 
flames,  LES  models  that  depend  on  scalar  gradients  are  susceptible  to  high  numerical  errors 
[99,  100],  which  would  partially  explain  the  over  prediction  here. 

Although  the  scalar  comparisons  upstream  (Fig.  48)  are  reasonably  good,  it  is  not 
known  if  the  minor  errors  seen  upstream  will  propagate  and  amplify  further  downstream. 
For  instance,  the  mixture  fraction  plots  show  an  over-prediction  of  roughly  15%  near  the 
centerline,  while  the  temperature  itself  is  nearly  the  same  as  experiments.  This  mismatch 
could  lead  to  fuel-rich  high  temperature  regions  downstream  that  will  promote  the  formation 
of  soot  particles.  Without  further  information  on  the  scalar  profiles  downstream,  it  is  not 
possible  to  deduce  the  role  of  these  errors  in  prediction  errors. 

While  comparing  experiments  to  simulations,  it  is  also  important  to  consider  the  exper¬ 
imental  errors.  Methane  flames  are  not  highly  sooting  flames,  with  the  peak  soot  volume 
fraction  in  the  parts  per  billion  levels.  Further,  the  structure  and  characteristics  of  the  soot 
particles  are  not  known  from  the  experiments.  Given  the  experimental  uncertainty  associ¬ 
ated  with  LH-based  measurement  of  nascent  low  volume  fraction  soot  signals,  it  is  likely 
that  these  errors  contribute  to  the  mismatch  as  well. 

4. 3. 4. 3  Subfilter  soot  and  scalar  statistics  From  the  particle  information,  subfilter 
statistics  and  PDFs  can  be  computed  to  better  understand  soot  evolution.  Figure  50  shows 
the  instantaneous  contours  of  subfilter  standard  deviation  of  mixture  fraction,  temperature, 
soot  volume  fraction,  and  number  density.  It  is  seen  that  the  subfilter  fluctuations  of  mixture 
fraction  are  confined  to  the  near-stoichiometric  region,  with  the  variations  essentially  reduced 
to  zero  outside  of  this  region.  While  temperature  fluctuations  follow  a  similar  pattern,  the 
high  standard  deviation  regions  are  seen  on  the  lean-side  of  the  flame,  outside  the  closed 
surface  formed  by  the  stoichiometric  contour.  It  should  be  noted  that  small  fluctuations 
in  mixture  fraction  could  lead  to  large  changes  in  temperature  close  to  the  stoichiometric 
region.  The  subfilter  fluctuations  of  soot  moments  are  also  large,  reaching  as  high  as  40%  of 
the  mean  values.  Similar  to  the  filtered  moments,  the  peaks  in  these  quantities  are  located 
in  the  relatively  high-temperature  fuel-rich  regions  inside  the  stoichiometric  surface. 

Figure  51  shows  the  subfilter  PDF  of  scalars  computed  at  different  axial  locations,  slightly 
off  the  centerline.  The  progress  variable  is  skewed  towards  the  values  greater  than  0.8, 
indicating  a  fully  burning  mixture.  The  mixture  fraction  PDF  shows  a  broader  spread,  with 
significant  probability  of  having  fuel- rich  pockets  (Z  >  0.075).  The  PAH  PDF  shows  a  bi- 
modal  behavior  with  a  large  fraction  of  the  particles  having  near-zero  PAH  concentrations 
and  another  group  of  particles  with  higher  PAH  values.  These  two  groups  of  particles  come 
from  the  fuel-lean  and  fuel-rich  sides,  respectively,  corresponding  to  the  broad  mixture  PDF 
observed.  The  PDF  of  soot  moments  exhibits  a  broad  and  nearly  uniform  PDF,  with  a  bias 
towards  zero  moment  corresponding  to  the  fuel-lean  mixture  fraction  values. 

4. 3. 4. 4  Conditional  statistics  The  evolution  of  the  soot  particles  in  mixture  fraction 
and  reaction  progress  variable  coordinates  is  illustrative  of  the  evolution  process.  Figure 
52  shows  the  one- dimensional  conditional  scatter  plots  of  PAH  mass  fraction,  soot  volume 
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Figure  50:  Instantaneous  contours  of  subfilter  standard  deviation  of  mixture  fraction,  tem¬ 
perature,  soot  volume  fraction  and  soot  number  density.  Iso- value  of  stoichiometric  mixture 
fraction  is  indicated  by  a  solid  black  line. 
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(d)  (e)  (f) 


Figure  51:  Time  averaged  marginal  PDFs  of  (a)  mixture  fraction,  (b)  reaction  progress 
variable,  (c)  temperature,  (d)  PAH  mass  fraction,  (e)  soot  mass  fraction,  and  (f)  soot 
number  density  plotted  at  r/d=1.5.  x/d=70:  (^—  ),  x/d=90:  (----  ) 


fraction  and  soot  number  density  at  x/d=70.  Here,  data  sets  at  several  instances,  spanning 
more  than  70  ms,  have  been  pooled  together.  All  three  quantities  show  a  peak  in  values 
near  a  mixture  fraction  of  0.165,  which  seems  to  be  most  likely  fuel  composition  to  contain 
soot.  It  is  also  seen  that  that  there  is  a  near  linear  increase  in  the  quantities  in  the  region 
of  0.075  <  Z  <  0.16,  but  profiles  have  a  non-linear  shape  in  richer  mixtures.  To  further 
understand  these  conditional  statistics,  Fig.  53  shows  the  instantaneous  conditional  scatter 
plot  at  two  different  locations  at  one  time  (same  as  that  in  Figs.  47  and  50).  It  is  seen 
that  the  when  the  filtered  soot  volume  fraction  is  lower  (<  20  ppb),  the  conditional  plots 
show  a  wider  spread  with  richer  regions  containing  significant  soot  volume  fractions.  Here, 
the  mixture  fraction  PDF  itself  is  broader,  with  a  wider  range  of  values  present.  However, 
when  the  filtered  soot  volume  fraction  is  higher  (x/d=90),  the  mixture  fraction  distribution 
appears  narrower,  and  the  conditional  scatter  plot  exhibits  a  linear  profile.  These  two  modes 
result  in  the  behavior  see  in  Fig.  52. 

Figure  54  shows  the  two-dimensional  scatter  plots  of  PAH  mass  fraction  and  soot  volume 
fraction  in  {Z,  C}  space.  The  plots  show  that  PAH  is  significantly  higher  in  a  small  region 
with  very  high  progress  variable  in  the  fuel-rich  region.  This  indicates  that  the  formation 
of  the  soot  precursors  is  confined  to  high  temperature  regions  on  the  fuel-rich  side  of  the 
stoichiometric  surface  (see  Fig.  47  ).  Soot  volume  fraction,  on  the  other  hand,  has  a  wider 
spread  with  large  soot  volume  fractions  seen  in  regions  of  lower  progress  variable  (and  lower 
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x/d=70 


Mixture  fraction  w  Mixture  fraction  Mixture  fraction 


Figure  52:  Conditional  scatter  plots  of  PAH  mass  fraction,  soot  volume  fraction  and  soot 
number  density  plotted  at  x/d=70.  Data  obtained  from  multiple  instantaneous  fields. 


temperature)  inside  the  fuel-rich  region.  This  indicates  that  the  soot  particles  travel  away 
from  the  region  of  high  PAH  concentration. 

It  is  also  important  to  note  that  these  conditional  scatter  plots  are  sensitive  to  the  subfilter 
mixing  model  used.  Since  soot  particles  do  not  diffuse,  the  assumption  of  equal  diffusivity 
will  have  an  impact  on  the  results  obtained  here.  In  particular,  soot  oxidation  in  the  near- 
stoichiometric  region  should  be  higher  than  that  see  in  this  simulation.  Since  the  IEM  model 
allows  the  notional  particles  to  mix  without  reacting,  allowing  traversal  in  composition  space 
across  flame  fronts,  the  effect  of  oxidation  is  under-predicted.  The  influence  of  differential 
diffusion  on  soot  evolution  needs  further  consideration 
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Mixture  fraction  Mixture  fraction 

Figure  53:  Instantaneous  variation  of  soot  volume  fraction  with  mixture  fraction  plotted  at 
x/d=70  (left)  and  x/d=90  (right). 


Figure  54:  Time  averaged  contours  of  PAH  mass  fraction  and  soot  volume  fraction  condi¬ 
tioned  on  mixture  fraction  and  reaction  progress  variable  plotted  at  x/d=70. 
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5  Conclusions 


5.1  Conclusions  on  Soot  Oxidation 

We  studied  the  thermodynamics  and  kinetics  of  aromatic  oxyradicals,  the  key  intermediates 
in  oxidation  of  soot. 

The  results  of  the  thermodynamics  study  show  that  thermodynamic  stability  varies  sub¬ 
stantially  at  different  chemisorbed  O  sites. 

Similarly,  the  results  of  the  kinetics  study  show  that  the  rate  of  oxyradical  decomposition 
varies  with  the  location  of  the  chemisorbed  O  site. 

The  principal  conclusion,  combining  the  thermodynamic  and  kinetic  results,  is  that  only 
the  “outer”  rings  are  able  to  undergo  oxidation,  similar  to  phenoxy  decomposition,  whereas 
the  inner  rings  resist  oxidation  in  combustion  environments.  This  finding  implies  that  for  an 
arbitrary  shaped  PAH,  oxidation  should  predominantly  remove  armchair  and  corner-zigzag 
sites,  leaving  resistant-to-oxidation  inner  zigzag  sites  essentially  intact.  Considering  that  the 
growth  of  both  armchair  and  zigzag  edges  proceed  at  effectively  the  same  rate,  we  expect  to 
find  proliferation  of  zigzag-edge  surfaces  on  soot  particles  formed  in  flame  environments. 

We  also  found  that  the  variability  of  both  thermodynamic  stability  and  decomposition 
rate  of  aromatic  oxyradicals  can  be  explained  and  correlated  with  substrate  aromaticity. 
These  correlations  could  lead  to  development  of  practical  rules  for  predicting  oxidation  rate 
coefficients  for  arbitrary  size  and  shape  aromatics,  thereby  supplementing  and  enhancing 
models  of  soot  oxidation. 

5.2  Conclusions  of  DNS  Soot  Study 

In  this  project,  we  carried  out  the  first  DNS  study  of  a  turbulent  sooting  flame  including 
finite  rate  chemistry  of  Polycyclic  Aromatic  Hydrocarbons  (PAH),  which  are  widely  accepted 
as  being  precursors  to  soot  particles. 

Our  results  on  the  formation  and  early  evolution  of  soot  highlight  several  features  having 
important  implications  in  turbulent  sooting  flames  research  and  modeling. 

•  The  concentration  of  PAH  is  a  function  of  local  mixture  fraction  and  scalar  dissipation 
rate.  Furthermore,  unsteady  mixing  appears  to  significantly  affect  PAH  concentration. 
This  behavior  explains  the  observed  sensitivity  of  soot  yield  to  strain  rates.  These 
findings  emphasize  that  turbulent  combustion  models  need  to  account  adequately  for 
the  dependence  of  PAH  concentrations  on  turbulent  mixing.  In  addition,  acetylene 
was  shown  to  be  less  sensitive  to  strain.  Semi-empirical  models  using  acetylene  to 
describe  soot  growth  are  not  expected  to  capture  the  correct  dependency  of  soot  yield 
on  mixing. 

•  Differential  diffusion  of  soot  due  to  negligible  molecular  diffusion  transport  contributes 
to  spreading  soot  in  mixture  fraction  space.  While  transport  processes  are  approxi¬ 
mately  symmetric  and  movement  towards  and  away  from  the  flame  is  equally  likely 
at  early  times,  the  distribution  of  soot  in  mixture  fraction  space  is  significantly  bi¬ 
ased.  The  highly  non-linear  growth  and  destruction  processes  affecting  soot  evolution 
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result  in  soot  being  preferentially  located  at  mixtures  richer  than  those  where  particles 
nucleate. 

•  Surface  reactions  are  most  active  next  to  the  flame  sheet  where  soot  is  oxidized  by 
hydroxyl  radicals  and  molecular  oxygen.  Soot  displaced  towards  the  flame  grows  by 
surface  reactions  and  is  oxidized  quickly.  In  regions  of  intense  growth,  the  displacement 
speed  is  usually  large,  leaving  little  time  for  growth. 

•  The  morphology  of  soot  particles  and  aggregates  is  largely  determined  by  their  motion 
in  mixture  fraction  space. 

5.3  Conclusions  on  LES  Soot  Model 

In  this  project,  an  integrated  Large  Eddy  Simulation  (LES)  model  for  soot  evolution  in 
turbulent  nonpremixed  flames  has  been  developed  and  validated  in  a  laboratory  scale  natural 
gas  jet  diffusion  flame.  The  model  components  include  an  accurate  yet  efficient  soot  model 
based  on  the  Hybrid  Method  of  Moments  (HMOM),  an  extended  Flamelet/Progress  Variable 
(RFPV)  gas-phase  combustion  model,  which  includes  heat  losses  due  to  gas-phase  radiation 
and  which  captures  the  detailed  chemical  kinetics  of  both  fuel  oxidation  and  soot  precursor 
formation,  and  a  model  to  account  for  the  slow  chemistry  of  PAH,  the  immediate  precursors 
of  soot.  Presumed  and  transported  PDF  approaches  for  subfilter  closure  of  the  combustion 
and  soot  models  were  considered. 

The  model  was  validated  in  Delft  Flame  III,  a  piloted  coflow  turbulent  nonpremixed  flame 
fueled  by  commercial  natural  gas.  Results  for  velocity,  mixture  fraction,  temperature,  and 
major  species  mass  fractions  in  the  upstream,  non-sooting  region  of  the  flame  showed  very 
good  agreement  with  experimental  measurements,  especially  compared  to  previous  RANS 
studies. 

Compared  to  the  experimental  measurements,  the  LES  provided  a  reasonable  prediction 
of  the  maximum  soot  volume  fraction.  This  quantity  was  found  to  be  sensitive  to  the  de¬ 
scription  of  small  scale  mixing  rates.  In  the  presumed  PDF  approach,  the  relevant  parameter 
is  the  subfilter  scalar  dissipation  rate.  Improved  agreement  was  obtained  by  decreasing  the 
model  constant  for  this  quantity,  which  increases  the  subfilter  scalar  variance.  Similarly, 
simulations  using  the  transported  PDF  were  found  to  be  sensitive  to  the  mixing  constant 
used  in  the  micromixing  model. 

However,  the  LES  predicted  the  formation  of  soot  much  closer  to  the  nozzle  than  in  the 
experiment.  The  primary  reason  for  this  discrepancy  was  attributed  to  significant  uncertainty 
in  the  chemistry  of  PAH.  The  LES  was  also  compared  to  experimental  measurements  of  the 
total  intermittency,  the  probability  of  soot  being  less  than  a  specified  threshold  at  a  given 
location  at  a  given  time.  While  quantitative  comparisons  are  not  meaningful  in  light  of  the 
disagreement  in  the  volume  fraction  profile,  the  LES  and  experimental  measurements  both 
showed  a  gradual  increase  in  the  intermittency  as  soot  is  oxidized.  This  increase  is  indicative 
of  a  decrease  in  the  frequency  of  pockets  of  large  soot  volume  fraction  in  the  soot  oxidation 
zone  rather  than  a  simple  decrease  in  the  soot  volume  fraction. 

While  this  work  is  major  step  forward  in  the  modeling  of  soot  evolution  with  LES,  two 
major  modeling  needs  that  have  a  significant  impact  on  the  results  presented  here  have  been 
identified.  First,  the  results  were  shown  to  be  quite  sensitive  to  the  model  coefficients  used  to 
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parameterize  small  scale  mixing  rates.  More  reliable  models  for  variance  and  dissipation  rate 
in  presumed  PDF  approaches  and  for  micromixing  in  transported  PDF  approaches  would 
lead  to  improvements  in  predictions  of  soot  evolution  in  turbulent  combustion.  Second, 
considerable  uncertainty  remains  in  the  kinetic  mechanism  for  PAH  formation  and  growth. 
Both  these  issues  will  be  addressed  in  future  research  efforts. 


77 


Bibliography 


[1]  P.  Stroomer.  Turbulence  and  OH  structures  in  flames.  PhD  thesis,  Delft  University  of 
Technology,  1995. 

[2]  P.  A.  Nooren,  M.  Versluis,  T.  H.  van  der  Meer,  R.  S.  Barlow,  and  J.  H.  Frank.  Raman- 
Rayleigh-LIF  measurements  of  temperature  and  species  concentrations  in  the  Delft 
piloted  turbulent  jet  diffusion  flame.  Appl.  Phys.  B ,  71:95-111,  2000. 

[3]  N.  H.  Qamar,  Z.  T.  Alwahabi,  Q.  N.  Chan,  G.  J.  Nathan,  D.  Roekaerts,  and  K.  D. 
King.  Soot  volume  fraction  in  a  piloted  turbulent  jet  non-premixed  flame  of  natural 
gas.  Combust.  Flame ,  156:1339-1347,  2009. 

[4]  M.  Frenklach.  Reaction  mechanism  of  soot  formation  in  flames.  Phys.  Chem.  Chem. 
Phys.,  4:2028-2037,  2002. 

[5]  H.  Bockhorn,  A.  D’Anna,  A.  F.  Sarofim,  and  H.  Wang.  Combustion  Generated  Fine 
Carbonaceous  Particles.  KIT  Scientific,  Karlsruhe,  Germany,  2009. 

[6]  H.  Wang.  Formation  of  nascent  soot  and  other  condensed-phase  materials  in  flames. 
Proc.  Combust.  Inst.,  33:41-67,  2011. 

[7]  M.  Frenklach.  A  Unifying  Picture  of  Gas  Phase  Formation  and  Growth  of  PAH,  Soot, 
Diamond  and  Graphite,  pages  259-273.  NASA  Conference  Publication  3061,  1990. 

[8]  M.  Frenklach  and  H.  Wang.  Detailed  modeling  of  soot  particle  nucleation  and  growth. 
Proc.  Combust.  Inst.,  23:1559-1566,  1991. 

[9]  M.  Frenklach.  On  surface  growth  mechanism  of  soot  particles.  Proc.  Combust.  Inst., 
26:2285-2293,  1996. 

[10]  R.  Whitesides  and  M.  Frenklach.  Detailed  kinetic  monte  carlo  simulations  of  graphene- 
edge  growth.  J.  Phys.  Chem.  A,  114:689-703,  2010. 

[11]  J.  Appel,  H.  Bockhorn,  and  M.  Frenklach.  Kinetic  modeling  of  soot  formation  with 
detailed  chemistry  and  physics,  laminar  premixed  flames  of  c2  hydrocarbons.  Combust. 
Flame,  121:122-136,  2000. 

[12]  K.  G.  Neoh,  J.  B.  Howard,  and  A.  F.  Sarofim.  page  261.  Plenum,  New  York,  1981. 

[13]  M.  S.  Celnik,  M.  Sander,  A.  Raj,  R.  H.  West,  and  M.  Kraft.  Modelling  soot  formation 
in  a  premixed  flame  using  an  aromatic-site  soot  model  and  an  improved  oxidation  rate. 
Proc.  Combust.  Inst.,  32:639-646,  2009. 


78 


[14]  B.  Zhao,  Z.  Yang,  M.V.  Johnston,  H.  Wang,  A.S.  Wexler,  M.  Balthasar,  and  M.  Kraft. 
Measurement  and  numerical  simulation  of  soot  particle  size  distribution  functions  in  a 
laminar  premixed  ethylene-oxygen-argon  flame.  Combustion  and  Flame ,  133:173-188, 
2003. 

[15]  B.  Zhao,  Z.  Yang,  Z.  Li,  M.V.  Johnston,  and  H.  Wang.  Particle  size  distribution  func¬ 
tion  of  incipient  soot  in  laminar  premixed  ethylene  flames:  effect  of  flame  temperature. 
Proceedings  of  the  Combustion  Institute ,  30:1441-1448,  2005. 

[16]  U.O.  Koylii,  GM  FaethT  L,  and  MG  Carvalho.  Fractal  and  projected  structure  prop¬ 
erties  of  soot  aggregates.  Combustion  and  Flame ,  100:621-633,  1995. 

[17]  D.L.  Marchisio  and  R.O.  Fox.  Solution  of  population  balance  equations  using  the  direct 
quadrature  method  of  moments.  Journal  of  Aerosol  Science,  36:43-73,  2005. 

[18]  M.  E.  Mueller,  G.  Blanquart,  and  H.  Pitsch.  A  joint  Volume-Surface  model  of  soot 
aggregation  with  the  method  of  moments.  Proc.  Combust.  Inst.,  32:785-792,  2009. 

[19]  G.  Blanquart  and  H.  Pitsch.  A  joint  Volume-Surface-Hydrogen  multi-variate  model 
for  soot  formation.  In  H.  Bockhorn,  A.  D’Anna,  A.  Sarofim,  and  H.  Wang,  editors, 
Combustion  Generated  Fine  Carbonaceous  Particles,  pages  439-466.  KIT  Scientific 
Publishing,  2009. 

[20]  M.  Frenklach  and  S.  J.  Harris.  Aerosol  dynamics  modeling  using  the  method  of  mo¬ 
ments.  J.  Coll.  Int.  Sci.,  118:252-261,  1987. 

[21]  M.  Frenklach.  Method  of  moments  with  interpolative  closure.  Chemical  Engineering 
Science,  57:2229-2239,  2002. 

[22]  R.  McGraw.  Description  of  aerosol  dynamics  by  the  quadrature  method  of  moments. 
Aerosol  Science  and  Technology,  27:255-265,  1997. 

[23]  M.  E.  Mueller,  G.  Blanquart,  and  H.  Pitsch.  Hybrid  Method  of  Moments  for  modeling 
soot  formation  and  growth.  Combust.  Flame,  156:1143-1155,  2009. 

[24]  F.  Bisetti,  G.  Blanquart,  M.  E.  Mueller,  and  H.  Pitsch.  On  the  formation  and  early 
evolution  of  soot  in  turbulent  nonpremixed  flames.  Combust.  Flame,  In  Press. 

[25]  M.  E.  Mueller,  G.  Blanquart,  and  H.  Pitsch.  Modeling  the  oxidation-induced  frag¬ 
mentation  of  soot  aggregates  in  laminar  flames.  Proc.  Combust.  Inst.,  33:667-674, 
2011. 

[26]  H.  Pitsch.  Large-eddy  simulation  of  turbulent  combustion.  Annual  Review  of  Fluid 
Mechanics,  38:453-482,  2006. 

[27]  H.  El-Asrag  and  S.  Menon.  Large  eddy  simulation  of  soot  formation  in  a  turbulent 
non-premixed  jet  flame.  Combust.  Flame,  156:385-395,  2009. 


79 


[28]  D.  O.  Lignell,  J.  H.  Chen,  P.  J.  Smith,  T.  Lu,  and  C.  K.  Law.  The  effect  of  flame 
structure  on  soot  formation  and  transport  in  turbulent  nonpremixed  flames  using  direct 
numerical  simulation.  Combust.  Flame ,  151:2-28,  2007. 

[29]  D.  O.  Lignell,  J.  H.  Chen,  and  P.  J.  Smith.  Three-dimensional  direct  numerical  simu¬ 
lation  of  soot  formation  and  transport  in  a  temporally  evolving  nonpremixed  ethylene 
jet  flame.  Combust.  Flame ,  155:316-333,  2008. 

[30]  M.  E.  Mueller  and  H.  Pitsch.  LES  subfilter  modeling  of  soot-turbulence  interactions. 
Phys.  Fluids,  Submitted. 

[31]  R.P.  Lindstedt  and  S.A.  Louloudi.  Joint-scalar  transported  pdf  modeling  of  soot  for¬ 
mation  and  oxidation.  Proceedings  of  the  Combustion  Institute,  30(1):775  -  783,  2005. 

[32]  Abhilash  J.  Chandy,  David  J.  Glaze,  and  Steven  H.  Frankel.  A  Hybrid  Large  Eddy 
Simulation/Filtered  Mass  Density  Function  for  the  Calculation  of  Strongly  Radiating 
Turbulent  Flames.  Journal  of  Heat  Transfer  -  Transactions  of  the  ASME,  131(5),  May 
2009. 

[33]  J.  P.  Perdew,  K.  Burke,  and  M.  Ernzerhof.  Generalized  gradient  approximation  made 
simple.  Phys.  Rev.  Lett.,  77:3865-3868,  1996. 

[34]  O.  F.  Sankey  and  D.  J.  Niklewski.  Abinitio  multicenter  tight-binding  model  for 
molecular-dynamics  simulations  and  other  applications  in  covalent  systems.  Phys. 
Rev.  B,  40:3979-3995,  1989. 

[35]  J.  M.  Soler,  E.  Artacho,  J.  D.  Gale,  A.  Garcia,  J.  Junquera,  P.  Ordejon,  and  D.  Sanchez- 
Portal.  The  siesta  method  for  ab  initio  order-n  materials  simulation.  J.  Phys.- 
Condensed  Matter,  14:2745-2779,  2002. 

[36]  J.  J.  P.  Stewart.  Optimization  of  parameters  for  semiempirical  methods  v:  Modification 
of  nddo  approximations  and  application  to  70  elements.  J.  Mol.  Modeling,  13:1173- 
1213,  2007. 

[37]  J  Kruszewski  and  T.  M.  Krygowski.  Definition  of  aromaticity  basing  on  harmonic 
oscillator  model.  Tetrahedron  Lett.,  36:3839,  1972. 

[38]  A.  D.  Becke.  Density-functional  thermochemistry  .3.  the  role  of  exact  exchange.  J. 
Chem.  Phys.,  98:5648-5652,  1993. 

[39]  R.  Krishnan,  J.  S.  Binkley,  R.  Seeger,  and  J.  A.  Pople.  Self-consistent  molecular-orbital 
methods  .20.  basis  set  for  correlated  wave- functions.  J.  Chem.  Phys.,  72:650-654,  1980. 

[40]  J.  R.  Barker.  Int.  J.  Chem.  Kinet.,  41:748-763,  2009. 

[41]  R.  G.  Gilbert  and  S.  C.  Smith.  Theory  of  Unimolecular  and  Recombination  Reactions. 
Blackwell-Scientific,  Oxford,  UK,  1990. 

[42]  H.  Wang  and  M.  Frenklach.  Combust.  Flame,  96:163-170,  1994. 


80 


[43]  D.  M.  Golden.  What  do  we  know  about  the  iconic  system  CH3+  CH3+  M  C2H6+ 
M?  Z.  Phys.  Chem.,  225:1-14,  2011. 

[44]  D.  L.  Marchisio  and  R.  O.  Fox.  Solution  of  population  balance  equations  using  the 
Direct  Quadrature  Method  of  Moments.  J.  Aerosol  Sci.,  36:43-73,  2005. 

[45]  C.  D.  Pierce  and  P.  Moin.  Progress  variable  approach  for  Large  Eddy  Simulation  of 
non-premixed  turbulent  combustion.  J.  Fluid  Mech.,  504:73-97,  2004. 

[46]  N.  Peters.  Laminar  diffusion  flamelet  models  in  non-premixed  turbulent  combustion. 
Prog.  Energy  Combust.  Sci.,  10:319-339,  1984. 

[47]  M.  Ihme  and  H.  Pitsch.  Modeling  of  radiation  and  nitric  oxide  formation  in  turbu¬ 
lent  nonpremixed  flames  using  a  flamelet /progress  variable  formulation.  Phys.  Fluids , 
20:055110,  2008. 

[48]  H.  Pitsch  and  N.  Peters.  A  consistent  flamelet  formulation  for  non-premixed  combus¬ 
tion  considering  differential  diffusion  effects.  Combust.  Flame,  114:26-40,  1998. 

[49]  R.  W.  Bilger.  The  structure  of  turbulent  nonpremixed  flames.  Proc.  Combust.  Inst., 
22:475-488,  1988. 

[50]  R.  S.  Barlow,  A.  N.  Karpetis,  J.  H.  Frank,  and  J.-Y.  Chen.  Scalar  profiles  and  NO 
formation  in  laminar  opposed-flow  partially  premixed  methane/air  flames.  Combust. 
Flame,  127:2102-2118,  2001. 

[51]  A.  W.  Cook  and  J.  J.  Riley.  A  subgrid  model  for  equilibrium  chemistry  in  turbulent 
flows.  Phys.  Fluids,  6:2868-2870,  1994. 

[52]  J.  Jimenez,  A.  Linan,  M.  M.  Rogers,  and  J.  F.  Higuera.  A  priori  testing  of  subgrid 
models  for  chemically  reacting  non-premixed  turbulent  shear  flows.  J.  Fluid  Mech., 
349:149-171,  1997. 

[53]  C.  Wall,  B.  J.  Boersma,  and  P.  Moin.  An  evaluation  of  the  assumed  beta  proba¬ 
bility  density  function  subgrid-scale  model  for  large  eddy  simulation  of  nonpremixed, 
turbulent  combustion  with  heat  release.  Phys.  Fluids ,  12:2522-2529,  2000. 

[54]  C.  D.  Pierce  and  P.  Moin.  A  dynamic  model  for  subgrid-scale  variance  and  dissipation 
rate  of  a  conserved  scalar.  Phys.  Fluids,  10:3041-3044,  1998. 

[55]  C.  M.  Kaul,  V.  Raman,  G.  Balarac,  and  H.  Pitsch.  Numerical  errors  in  the  computation 
of  subfilter  scalar  variance  in  large  eddy  simulations.  Phys.  Fluids,  21:055102,  2009. 

[56]  C.  M.  Kaul  and  V.  Raman.  A  posteriori  analysis  of  numerical  errors  in  subfilter  scalar 
variance  modeling  for  large  eddy  simulation.  Phys.  Fluids,  23:035102,  2011. 

[57]  M.  Ihme  and  H.  Pitsch.  Prediction  of  extinction  and  reignition  in  nonpremixed  turbu¬ 
lent  flames  using  a  flamelet /progress  variable  model.  2.  Application  in  LES  of  Sandia 
flames  D  and  E.  Combust.  Flame,  155:90-107,  2008. 


81 


[58]  S.  James,  J.  Zhu,  and  M.  S.  Anand.  Lagrangian  PDF  transport  method  for  simulations 
of  axisymmetric  turbulent  reacting  flows.  In  43rd  AIAA  Aerospace  Sciences  Meeting 
and  Exhibit ,  number  156,  2005. 

[59]  S.  B.  Pope.  A  Monte  Carlo  method  for  the  PDF  equations  of  turbulent  reactive  flow. 
Combustion  Science  and  Technology ,  25:159-174,  1981. 

[60]  P.  J.  Colucci,  F.  A.  Jaberi,  and  P.  Givi.  Filtered  density  function  for  large  eddy 
simulation  of  turbulent  reacting  flows.  Physics  of  Fluids,  10(2):499-515,  1998. 

[61]  V.  Raman,  H.  Pitsch,  and  R.  O.  Fox.  A  consistent  hybrid  LES-FDF  scheme  for  the 
simulation  of  turbulent  reactive  flows.  Combustion  and  Flame,  143(1-2) :56-78,  2005. 

[62]  H.  Pitsch,  E.  Riesmeier,  and  N.  Peters.  Unsteady  flamelet  modeling  of  soot  formation 
in  turbulent  diffusion  flames.  Combustion  science  and  technology ,  158:389-406,  2000. 

[63]  Y.  Liu  and  R.  O.  Fox.  Cfd  predictions  for  chemical  processing  in  a  confined  impinging- 
jets  reactor.  AIChE  Journal,  52:731,  2006. 

[64]  R.  O.  Fox.  Computational  Models  for  Turbulent  Reacting  Flows.  Cambridge  University 
Press,  Cambridge,  2003. 

[65]  F.  A.  Jaberi,  P.  J.  Colucci,  S.  James,  P.  Givi,  and  S.  B.  Pope.  Filtered  mass  den¬ 
sity  function  for  large-eddy  simulation  of  turbulent  reacting  flows.  Journal  of  Fluid 
Mechanics,  401:85-121,  1999. 

[66]  V.  Raman  and  H.  Pitsch.  A  consistent  LES/filtered-density  function  formulation  for 
the  simulation  of  turbulent  flames  with  detailed  chemistry.  Proceedings  Proceedings  of 
the  Combustion  Institute,  31:1711-1719,  2006. 

[67]  H.  Pitsch.  FlameMaster.  A  C++  computer  program  for  0D  combustion  and  ID  laminar 
flame  calculations. 

[68]  G.  Blanquart,  P.  Pepiot-Desjardins,  and  H.  Pitsch.  Chemical  mechanism  for  high 
temperature  combustion  of  engine  relevant  fuels  with  emphasis  on  soot  precursors. 
Combust.  Flame,  156:588-607,  2009. 

[69]  K.  Narayanaswamy,  G.  Blanquart,  and  H.  Pitsch.  A  consistent  chemical  mechanism 
for  oxidation  of  substituted  aromatic  species.  Combust.  Flame,  157:1879-1898,  2010. 

[70]  C.  S.  McEnally  and  L.  D.  Pfefferle.  An  experimental  study  in  non-premixed  flames 
of  hydrocarbon  growth  processes  that  involve  five-membered  carbon  rings.  Combust. 
Sci.  Tech.,  131:323-344,  1998. 

[71]  A.  Laskin  and  A.  Lifshitz.  Thermal  decomposition  of  indene:  Experimental  results 
and  kinetic  modeling.  Proc.  Combust.  Inst.,  27:313-320,  1998. 

[72]  M.  Frenklach,  C.  A.  Schuetz,  and  J.  Ping.  Migration  mechanism  of  aromatic-edge 
growth.  Proc.  Combust.  Inst.,  30:1389-1396,  2005. 


82 


[73]  R.  Whitesides,  D.  Domin,  R.  Salomon-Ferrer,  Jr.  Lester,  W.  A.,  and  M.  Frenklach. 
Embedded-ring  migration  on  graphene  zigzag  edge.  Proc.  Combust.  Inst.,  32(1):577- 
583,  2009. 

[74]  X.Q.  You,  R.  Whitesides,  D.  Zubarev,  Jr.  Lester  W.  A.,  and  M.  Frenklach.  Bay¬ 
capping  reactions:  Kinetics  and  influence  on  graphene-edge  growth.  Proc.  Combust. 
Inst.,  33:685-692,  2011. 

[75]  D.  Y.  Zubarev,  N.  Robertson,  D.  Domin,  J.  McClean,  J.  Wang,  Jr.  Lester,  W.  A., 
R.  Whitesides,  X.  You,  and  M.  Frenklach.  Local  electronic  structure  and  stability  of 
pentacene  oxyradicals.  J.  Phys.  Chem.  C,  114:5429-5437,  2010. 

[76]  D.  Y.  Zubarev,  X.  You,  J.  McClean,  Jr.  Lester,  W.  A.,  and  M.  Frenklach.  Patterns  of 
local  aromaticity  in  graphene  oxyradicals.  J.  Mater.  Chem.,  21:3404-3409,  2011. 

[77]  D.  Y.  Zubarev,  X.  You,  J.  McClean,  Jr.  Lester,  W.  A.,  and  M.  Frenklach.  Delocalization 
effects  in  pristine  and  oxidized  graphene  substrates.  In  P.  E.  Hoggan,  E.  J.  Brondas, 
J.  Maruani,  P.  Piecuch,  ,  and  G.  Delgado-Barrio,  editors,  Advances  in  the  Theory 
of  Quantum  Systems  in  Chemistry  and  Physics,  pages  553-569.  Springer,  Dordrecht, 
2012. 

[78]  X.  You,  D.  Yu.  Zubarev,  Jr.  Lester,  W.  A.,  and  M.  Frenklach.  Thermal  decomposition 
of  pentacene  oxyradicals.  J.  Phys.  Chem.  A,  115:14184-14190,  2011. 

[79]  L.  R.  Radovic.  Active  sites  in  graphene  and  the  mechanism  of  C02  formation  in  carbon 
oxidation.  J.  Am.  Chem.  Soc.,  131:17166-17175,  2009. 

[80]  N.  Kunioshi,  M.  Touda,  and  S.  Fukutani.  Computational  study  on  the  formation  of 
five-membered  rings  in  PAH  through  reaction  with  02.  Combust.  Flame,  128(3)  :292- 
300,  2002. 

[81]  C.  Y.  Lin  and  M.  C.  Lin.  Thermal  decomposition  of  methyl  phenyl  ether  in  shock 
waves:  the  kinetics  of  phenoxy  radical  reactions.  J.  Phys.  Chem.,  90(3):425-431,  1986. 

[82]  P.  Frank,  J.  Herzler,  Th  Just,  and  C.  Wahl.  High-temperature  reactions  of  phenyl 
oxidation.  Symp.  Int.  Combust.,  25:833-840,  1994. 

[83]  Norbert  Peters.  Turbulent  Combustion.  Cambridge  University  Press,  2000. 

[84]  H.  Pitsch,  E.  Riesmeier,  and  N.  Peters.  Unsteady  flamelet  modeling  of  soot  formation 
in  turbulent  diffusion  flames.  Combust.  Sci.  Tech.,  158:389-406,  2000. 

[85]  T.  W.  J.  Peeters,  P.  P.  J.  Stroomer,  J.  E.  de  Vries,  D.  J.  E.  M.  Roekaerts,  and  C.J. 
Hoogendoorn.  Comparative  experimental  and  numerical  investigation  of  a  piloted 
turbulent  natural-gas  diffusion  flame.  Proc.  Combust.  Inst.,  25:1241-1248,  1994. 

[86]  O.  Desjardins,  G.  Blanquart,  G.  Balarac,  and  H.  Pitsch.  High  order  conservative  finite 
difference  scheme  for  variable  density  low  mach  number  turbulent  flows.  J.  Comp. 
Phys.,  227:7125-7159,  2008. 


83 


[87]  M.  Hermann,  G.  Blanquart,  and  V.  Raman.  A  bounded  QUICK  scheme  for  preserving 
scalar  bounds  in  large-eddy  simulations.  AIAA  J.,  44:2879-2886,  2006. 

[88]  M.  Germano,  U.  Piomelli,  P.  Moin,  and  W.  H.  Cabot.  A  dynamic  subgrid-scale  eddy 
viscosity  model.  Phys.  Fluids  A,  3:1760-1765,  1991. 

[89]  P.  Moin,  K.  Squires,  W.  Cabot,  and  S.  Lee.  A  dynamic  subgrid-scale  model  for  com¬ 
pressible  turbulence  and  scalar  transport.  Phys.  Fluids  A,  3:2746-2757,  1991. 

[90]  D.  K.  Lilly.  A  proposed  modification  of  the  Germano  subgrid-scale  closure  method. 
Phys.  Fluids  A,  4:633-635,  1992. 

[91]  C.  Meneveau,  T.  S.  Lund,  and  W.  H.  Cabot.  A  Lagrangian  dynamic  subgrid-scale 
model  of  turbulence.  J.  Fluid  Mech.,  319:353-385,  1996. 

[92]  J.  Reveillon  and  L.  Vervisch.  Response  of  the  dynamic  LES  model  to  heat  release 
induced  effects.  Phys.  Fluids ,  8:2248-2250,  1996. 

[93]  B.  Merci,  B.  Naud,  and  D.  Roekaerts.  Flow  and  mixing  fields  for  transported  scalar 
PDF  simulations  of  a  piloted  jet  diffusion  flame  (’Delft  Flame  IIP).  Flow  Turb.  Com¬ 
bust.,  74:239-272,  2005. 

[94]  B.  Merci,  D.  Roekaerts,  and  B.  Naud.  Study  of  the  the  performance  of  three  micromix¬ 
ing  models  in  transported  scalar  PDF  simulations  of  a  piloted  jet  diffusion  flame  (’’Delft 
Flame  III”).  Combust.  Flame ,  144:476-493,  2006. 

[95]  B.  Merci,  B.  Naud,  and  D.  Roekaerts.  Interaction  between  chemistry  and  micro- mixing 
modeling  in  transported  PDF  simulations  of  turbulent  non-premixed  flames.  Combust. 
Set.  Tech. ,  179:153-172,  2007. 

[96]  A.  Habibi,  B.  Merci,  and  D.  Roekaerts.  Turbulence  radiation  interaction  in  Reynolds- 
averaged  Navier-Stokes  simulations  of  nonpremixed  piloted  turbulent  laboratory-scale 
flames.  Combust.  Flame ,  151:303-320,  2007. 

[97]  E.  H.  van  Veen  and  D.  Roekaerts.  On  the  accuracy  of  temperature  measurements  in 
turbulent  jet  diffusion  flames  by  Coherent  Anti-Stokes-Raman  spectroscopy.  Combust. 
Sci.  Tech.,  175:1893-1914,  2003. 

[98]  M.  E.  Mueller  and  H.  Pitsch.  LES  model  for  sooting  turbulent  nonpremixed  flames. 
Combust.  Flame,  In  Press.,  2012. 

[99]  C.  M.  Kaul,  V.  Raman,  G.  Balarac,  and  H.  Pitsch.  Numerical  errors  in  the  computation 
of  subfilter  scalar  variance  in  large  eddy  simulations.  Physics  of  Fluids,  21:055102,  2009. 

[100]  C.  M.  Kaul  and  V.  Raman.  A  posteriori  analysis  of  numerical  errors  in  subfilter  scalar 
variance  modeling  for  large  eddy  simulation.  Physics  of  Fluids,  23:035102,  2011. 


84 


List  of  Publications 

2009 

G.  Blanquart,  M.  E.  Mueller,  and  H.  Pitsch,  Modeling  and  Analysis  of  the  Effects  of  Tem¬ 
perature  on  Soot  Formation  using  a  Joint  Volume-Surface-Hydrogen  Model,  Comb. 
Flame,  156:1614-1626,  2009. 


G.  Blanquart,  P.  Pepiot-Desjardins,  and  H.  Pitsch,  Chemical  Mechanism  for  High  Temper¬ 
ature  Combustion  of  Engine  Relevant  Fuels  with  Emphasis  on  Soot  Precursors,  Comb. 
Flame,  156:588-607,  2009. 


M.  E.  Mueller,  G.  Blanquart,  and  H.  Pitsch,  Hybrid  Method  of  Moments  for  Modeling  Soot 
Formation  and  Growth,  Comb.  Flame,  156:1143-1155,  2009. 


M.  E.  Mueller,  G.  Blanquart,  and  H.  Pitsch,  A  Joint  Volume-Surface  Model  of  Soot  Aggre¬ 
gation  with  the  Method  of  Moments,  Proc.  Comb.  Inst.,  32:785-792,  2009. 


V.  Raman,  M.  E.  Mueller,  G.  Blanquart,  and  H.  Pitsch,  LES/PDF  modeling  of  soot  evo¬ 
lution  in  turbulent  flames,  62nd  Annual  Meeting  of  the  American  Physical  Society 
Division  of  Fluid  Dynamics,  Minneapolis,  MN,  November  22-24,  2009. 


2010 

J.  Wang,  D.  Domin,  B.  Austin,  D.  Zubarev,  J.  McClean,  M.  Frenklach,  T.  Cui,  and  W.  A. 
Lester,  Jr.,  A  Diffusion  Monte  Carlo  Study  of  the  O-H  Bond  Dissociation  of  Phenol, 
J.  Phys.  Chem.  A  114:9832-9835,  2010. 

D.  Y.  Zubarev,  N.  Robertson,  D.  Domin,  J.  McClean,  J.  Wang,  Jr.  Lester,  W.  A.,  R. 
Whitesides,  X.  You,  and  M.  Frenklach,  Local  electronic  structure  and  stability  of 
pentacene  oxyradicals,  J.  Phys.  Chem.  C,  114:5429-5437,  2010. 

M.  E.  Mueller  and  H.  Pitsch,  LES  Subfilter  Modeling  of  Soot-Turbulence  Interactions,  63rd 
Annual  Meeting  of  the  American  Physical  Society  Division  of  Fluid  Dynamics,  Long 
Beach,  CA,  November  21-23,  2010. 


2011 

F.  Bisetti,  G.  Blanquart,  M.  E.  Mueller,  and  H.  Pitsch,  On  the  formation  and  early  evolu¬ 
tion  of  soot  in  turbulent  nonpremixed  flames.  Combust.  Flame,  In  Press,  2011. 


M.  E.  Mueller,  G.  Blanquart,  and  H.  Pitsch,  Modeling  the  Oxidation-induced  Fragmenta¬ 
tion  of  Soot  Aggregates  in  Laminar  Flames.  Proc.  Combust.  Inst.,  33:667-674,  2011. 


85 


X.Q.  You,  R.  Whitesides,  D.  Zubarev,  Jr.  Lester  W.  A.,  and  M.  Frenklach,  Bay-capping 
reactions:  Kinetics  and  influence  on  graphene-edge  growth,  Proc.  Combust.  Inst., 
33:685-692,  2011. 

D.  Y.  Zubarev,  X.  You,  J.  McClean,  Jr.  Lester,  W.  A.,  and  M.  Frenklach,  Patterns  of  local 
aromaticity  in  graphene  oxyradicals,  J.  Mater.  Chem.,  21:3404-3409,  2011. 


2012 

D.  Y.  Zubarev,  X.  You,  J.  McClean,  Jr.  Lester,  W.  A.,  and  M.  Frenklach,  Delocalization 
effects  in  pristine  and  oxidized  graphene  substrates.  In  Advances  in  the  Theory  of 
Quantum  Systems  in  Chemistry  and  Physics  (P.  E.  Hoggan,  E.  J.  Brbndas,  J.  Maruani, 
P.  Piecuch,  and  G.  Delgado- Barrio,  Eds.),  pages  553-569.  Springer,  Dordrecht,  2012. 

P.  Donde,  V.  Raman,  M.  E.  Mueller,  and  H.  Pitsch,  LES/PDF  modeling  of  soot  evolution 
in  turbulent  flames,  Submitted  to  the  Proceedings  of  the  Combustion  Institute  as  part 
of  the  34th  Combustion  Symposium. 

M.  E.  Mueller  and  H.  Pitsch,  LES  model  for  sooting  turbulent  nonpremixed  flames,  Com¬ 
bustion  and  Flame,  In  Press,  2012. 


86 


